*! Ver 1.71 bug with weights
* Ver 1.7 adds version for easier update
* Ver 1.68 Bug with drimp
* Ver 1.67 Bug with IPW 
* Ver 1.66 Bug with binit, skip fixed
* Ver 1.65 Bug with binit. Fixed
* Ver 1.64 added binit for faster csdid (if slighly)
* Ver 1.63 added Dryrun
* Ver 1.63 Bug with option all and RC estimator
* Ver 1.62 Added level
* Ver 1.61 Change sample for drimp
* Ver 1.6 Change in check for 2x2 data. And updated site!
* added correction to teffects for ALL
* Ver 1.5   New output, with panel GMM estimators. Also WB standard errors with cluster
* Ver 1.38  Adding Cluster Stantandard errors
* version 1.37 2jun2021 Add extra messages of 2x2 balance. Checks that you indeed have panel data
* version 1.36 17may2021 Changes with EP display options
** Goal. Make the estimator modular. That way other options can be added
* v1.35 DRDID for Stata 
* Added by AN and adapted by FRA
* Added Version control (v14), Variables are no longer left after DRDID. But option stub is added 
* for later, as CSDID will use those
* v1.31 DRDID for Stata by FRA. Adding program for Bootstrap multiplier
* v1.3 DRDID for Stata by FRA Changing ATT RIF creation
* v1.2 DRDID for Stata by FRA All Estimators are ready For panel and RC
* Next help!
* v1.0 DRDID for Stata by FRA All Estimators but IPWRA have are available for panel
* Need to add weights. 
* v0.8 DRDID for Stata by FRA Other estimators are available. Onlyone missing iwp panel
* v0.7 DRDID for Stata by FRA IPW estimator for panel (Asjad original Rep)
* v0.5 DRDID for Stata by FRA Incorporates RC1 and RC2 estimators
* v0.2 DRDID for Stata by FRA Fixes typo with tag
* v0.2 DRDID for Stata by FRA Allows for Factor notation
* v0.1 DRDID for Stata by FRA Typo with ID TIME
/* !! epg: added method gmm . Also added parsing routing _Vce_Parse
      gmm and vce(...)                     
*/ 
/* !! epg: what to do with bootstrap csboot(csboot_opts)
				csboot_opts -- reps() 
				            -- rseed()
							-- wbtype()
							-- ...
							Perfecto ya todo fue integrado
  !! Crear la tabla Funciona
  !! Necesitamos regresar e(V) o solamente e(b)?
     Creo que ya esta. 
  !! Deberia cambiar el t estadistico?
	 
*/
/*
csdid loop over 
*/

program define drdid, eclass byable(onecall)
        version 14
        if _by() {
                local BY `"by `_byvars'`_byrc0':"'
        }

        `BY' _vce_parserun drdid, noeqlist jkopts(eclass): `0'
        
        if "`s(exit)'" != "" {
                ereturn local cmdline `"drdid `0'"'
                exit
        }
		
		syntax [anything(everything)] [iw pw aw], [* version]
		/**Version**/
		if   "`version'"!="" {
			display "version: 1.71"
			addr scalar version = 1.71
			exit
		}
		
        if replay() {
                if `"`e(cmd)'"' != "drdid" { 
                        error 301
                }
                else if _by() { 
                        error 190 
                }
                else {
                        Display `0'
                }
                exit
        }
		if runiform()<.001 {
			easter_egg
		}
        `vv' ///
        `BY' drdid_wh `0'
        ereturn local cmdline `"drdid `0'"'
		
end

mata
void is_2x2(string scalar time, treat, touse, isok){ 
    real matrix dta
 	
	dta=st_data(.,time,touse),st_data(.,treat,touse)
	dta=uniqrows(dta)
	if ((rows(dta)==4) & (rows(uniqrows(dta[,1]))==2) & (rows(uniqrows(dta[,2]))==2)) st_numscalar(isok,1)
	else st_numscalar(isok,0)
}
	

void is_balp(string scalar ivar, touse, balp){ 
    real matrix dta
 	dta=st_data(.,ivar,touse)
	st_numscalar(balp,mean(uniqrows(dta,2)[,2]))
} 
end

program define drdid_wh, eclass sortpreserve byable(recall)
	syntax varlist(fv numeric) [if] [in] [iw],			///
							[Ivar(varname)] 			///
							Time(varname) 				///
							TReatment(varname) 			///
							[noisily 					///
							drimp 						///
							dripw 						///
							reg 						///
							stdipw 						///
							ipw 						///
							ipwra 						///
							all  						///
							rc1 						///
							WBOOT(string) 				///
							WBOOT1						///
							*reps(int 999) 				///
							*wbtype(int 1)  			/// Hidden option
							rseed(str)					/// set seed
							Level(int 95)				/// CI level
							stub(name) replace 			/// to avoid overwritting
							cluster(varname)			/// For Cluster
							vce(string)					///
							gmm							///
							pscore(string)				///
							csdid						///
							binit(string)				///
							dryrun						///
							*							///
							]  
    * so it returns Nothing								
	if "`dryrun'"!="" error 1111
	
	_get_diopts diopts other, `options' 
	quietly capture Display, `diopts' `other' level(`level')	
	if _rc==198 {
		Display, `diopts' `other'  level(`level')     
	}


	marksample touse
	markout `touse' `ivar' `time' `treatment'  `cluster'
	
	_Vce_Parse if `touse',  `gmm' `wboot1' wboot(`wboot') vce(`vce')
** Move this into VCE_parse
	*if "`binit'"!="" 	local  binit `binit', skip
	local semethod "`r(semethod)'"
	
	if ("`semethod'"=="wildboot") {
	    local wboot "wboot"
		local reps = r(reps)
		local wbtype = r(wbtype)
	}
	
	 ** Verify 2x2 data
	tempname isok
	mata:is_2x2("`time'","`treatment'","`touse'","`isok'")
	if scalar(`isok')==0 {
	    display in red "You do not have a 2x2 design."
		error 555
	}	
			
			
	if "`cluster'"==""  local cluster "`r(cluster)'"
		
	if "`cluster'"!=""  {
		tempvar clvar
		qui:egen double `clvar' = group(`cluster') if `touse'
		local ocluster `cluster'
		local cluster `clvar'
	}
	
	if "`rseed'"=="" 	local rseed    "`r(seed)'"
	
	capture:set seed `rseed'
	
 	**# Verifies if data is panel data when "ivar" is declared
	if "`ivar'"!="" {
		capture xtset
		if _rc!=0 {
			qui:xtset `ivar' `time'
			qui:xtset ,  clear
		}
		if "`cluster'"!="" {
			_xtreg_chk_cl2 `cluster' `ivar'
		}
		
		**# verify ALL data is locally balanced
		tempname balp
 		*mata:is_balp("`ivar'", "`touse'", "`balp'")
		qui:bysort `touse' `ivar':gen `balp'=_N if `touse'
		sum `balp', meanonly
		if r(mean)<2 {
			display in red "{p}Some panel units are observed only once.{p_end}" _n ///
						   "{p}Those observations will be excluded from the sample. " ///
						   "If you want to keep them, do not use `ivar' as the panel id -ivar-, " ///
						   "which will request repeated crossection estimators {p_end}"
			quietly replace `touse' = 0 if `balp'==1
		}
	}
	
**# Verifies If variable exists. For CSDID
		
	cap unab allnew: `stub'att* 
	if "`stub'"!="" & "`replace'"=="" & "`allnew'"!="" {
		foreach v of var `allnew' {
			conf new var `v'
		}
	}
	
	*** Add Char to variables. They will be ours. Perhaps for CSDID
**# Verifies REPS is a number
	

	** Which option 
	if "`drimp'`dripw'`reg'`stdipw'`ipw'`ipwra'`all'"=="" {
	    display "No estimator selected. Using default {bf:drimp}"
		local drimp drimp
	}
	else {
	    if `:word count `drimp' `dripw' `reg' `stdipw' `ipw' `ipwra' `all' '!=1 {
		    display "Only one option allowed, more than 1 selected."
			error 1
		}
	}

	_S__est, `drimp' `dripw' `reg' `stdipw' `ipw' `ipwra' `all'
	local estimator "`s(estimator)'"

	** First determine outcome and xvars
	gettoken y xvar:varlist
	** Sanity Checks for Time. Only 2 values
	** just in case for xvar not be empty
	*local xvar `xvar'
	tempvar vals vals2 vals3
 
	
	tempvar tmt
	qui:egen byte `tmt'=group(`time') if `touse'
	qui:replace `tmt'=`tmt'-1	    
	tempvar trt
	qui:egen byte `trt'=group(`treatment') if `touse'
	qui:replace `trt'=`trt'-1
	
	** Sanity Check Weights
	** Weights
	if "`exp'"=="" {
	    tempname wgt
		gen byte `wgt'=1
	}
	else {
		tempvar wgt
		qui:gen double `wgt'`exp'
		qui:sum `wgt' if `touse' & `tmt'==0, meanonly
		qui:replace `wgt'=`wgt'/r(mean)
	}
	
	
	**# Here we collect all options
	** This are all the options send to DRDID 
	
	local 01 touse(`touse') tmt(`tmt') trt(`trt') y(`y') 	///
			 xvar(`xvar') `isily' ivar(`ivar') 	///
			 weight(`wgt') stub(`stub') ///
			  treatvar(`treatment') `rc1' cluster(`cluster') ///
			 `wboot' reps(`reps')  wbtype(`wbtype') level(`level')  binit(`binit')
			 *seed(`seed') 
		
	** Default will be IPT 
 	if "`estimator'"!="all" {
		if ("`semethod'"!="gmm") {
			drdid_`estimator', `01' `diopts' 
			ereturn local semethod `semethod'
			ereturn local clustvar `ocluster'
			ereturn local seed     `rseed'
			if "`ivar'"!="" ereturn local datatype "panel"
			else            ereturn local datatype "rcs"
			Display, bmatrix(e(b)) vmatrix(e(V)) `diopts' level(`level')
			exit 
		}
		else {
			tempvar touse2
			tempname b V 
			
			if ("`ivar'"=="") {
				local repeated repeated
				if ("`estimator'"=="sipwra") {
					display as error ///
					"{bf:ipwra} not allowed with {bf:gmm} and repeated"	///
					" cross-sectional data"	
					exit 198
				}
			}
			
			if "`weight'" != "" {
				local wgtgmm [`weight' `exp']
			}
			quietly generate byte `touse2' = . if `touse'

			_het_did_gmm `y' `xvar' if `touse' `wgtgmm', 				///
						  estimator(`estimator') groupvar(`ivar')		///
						  psvars(`xvar') treatvar(`trt')				///
						  timevar(`time') vce(`vce') touse2(`touse2')	///
						  t0(`tmt') pscore(`pscore') `repeated'			///
						  treatname(`treatment')

			matrix `b' = r(b)
			matrix `V' = r(V)
			local vce "`r(vce)'"
			local vcetype "`r(vcetype)'"
			local N = e(N)
			if ("`vce'"=="cluster") {
				local N_clust = r(N_clust)
				local clustvar "`r(clustvar)'"
			}
			ereturn post `b' `V', buildfvinfo esample(`touse2') obs(`N')
			ereturn local cmd drdid
			ereturn local method  `drimp'`dripw'`reg'`stdipw'`aipw'`ipwra'`all'
			ereturn local semethod `semethod'
			
			if ("`vce'"=="cluster") {
				ereturn scalar N_clust = `N_clust'
				ereturn local clustvar "`clustvar'"
			}
			ereturn local vce "`vce'"
			ereturn local vcetype "`vcetype'"
			ereturn local policy "`treatment'"
			ereturn hidden local method2 ///
				`drimp'`dripw'`reg'`stdipw'`aipw'`ipwra'`all'	
			Display, bmatrix(e(b)) vmatrix(e(V)) `diopts' 
			exit 
		}
	
	
	}
	 
	if "`estimator'"=="all" {
	** DR
	tempvar ttouse
	    qui:clonevar `ttouse'=`touse'
		tempname bb VV
		qui:drdid_dripw , `01'
 		matrix `bb' = nullmat(`bb')\e(b)
		matrix `VV' = nullmat(`VV')\e(V)
		local clname ATET:dripw	
		**capture drop `stub'att_dripw
		*ren `stub'att  `stub'att_dripw
		if "`ivar'"=="" {
			qui:capture gen byte `touse'=`ttouse'
			qui:drdid_dripw , `01' rc1 
			matrix `bb' = nullmat(`bb')\e(b)
			matrix `VV' = nullmat(`VV')\e(V)
			local clname `clname' ATET:dripw_rc1
			**capture drop `stub'att_dripwrc
			*clonevar `stub'att_dripwrc=`stub'att
		}
	** DRIMP	
		qui:capture gen byte `touse'=`ttouse'
 	    qui:drdid_imp , `01'
		matrix `bb' = nullmat(`bb')\e(b)
		matrix `VV' = nullmat(`VV')\e(V)
		local clname `clname' ATET:drimp
		**capture drop `stub'att_drimp
		*ren `stub'att  `stub'att_drimp
		
		if "`ivar'"=="" {
			qui:capture gen byte `touse'=`ttouse'
		    qui:drdid_imp , `01' rc1 
			matrix `bb' = nullmat(`bb')\e(b)
			matrix `VV' = nullmat(`VV')\e(V)
			local clname `clname' ATET:drimp_rc1
			**capture drop `stub'att_drimprc
			*ren `stub'att  `stub'att_drimprc
		}
	** REG
		qui:capture gen byte `touse'=`ttouse'
		qui:drdid_reg , `01'
		matrix `bb' = nullmat(`bb')\e(b)
		matrix `VV' = nullmat(`VV')\e(V)
		local clname `clname' ATET:reg
		**capture drop `stub'att_reg
		*ren `stub'att  `stub'att_reg
	** TRAD_IPW	
		qui:capture gen byte `touse'=`ttouse'
	    qui:drdid_aipw , `01'
		matrix `bb' = nullmat(`bb')\e(b)
		matrix `VV' = nullmat(`VV')\e(V)
		local clname `clname' ATET:ipw
		**capture drop `stub'att_ipw
		*ren `stub'att  `stub'att_ipw
	** STD IPW	
		qui:capture gen byte `touse'=`ttouse'
		qui:drdid_stdipw , `01'
		matrix `bb' = nullmat(`bb')\e(b)
		matrix `VV' = nullmat(`VV')\e(V)
		**capture drop `stub'att_stdipw
		*ren `stub'att  `stub'att_stdipw
		local clname `clname' ATET:stdipw
		if "`ivar'"!="" {
			qui:capture gen byte `touse'=`ttouse'
			qui:drdid_sipwra , `01'
			matrix `bb' = nullmat(`bb')\e(b)
			matrix `VV' = nullmat(`VV')\e(V)
			local clname `clname' ATET:sipwra
		}
		matrix `bb'=`bb''
		matrix colname `bb' = `clname'
		matrix `VV'=diag(`VV')
		matrix colname `VV' = `clname'		
		matrix rowname `VV' = `clname'
		ereturn post `bb' `VV'
		*ereturn display
	}
	
	ereturn local cmd drdid
	ereturn local policy "`treatment'"
	ereturn local method         `drimp'`dripw'`reg'`stdipw'`aipw'`ipwra'`all'
	ereturn hidden local method2 `drimp'`dripw'`reg'`stdipw'`aipw'`ipwra'`all'	

    Display, bmatrix(e(b)) vmatrix(e(V)) `diopts' level(`level')
end

program define _Vce_Parse, rclass
	syntax [anything] [if] [in], [gmm WBOOT1 vce(string) WBOOT(string)]
	
	marksample touse 
	
	local semethod standard
	if (("`wboot'"!=""|"`wboot1'"!="") & "`vce'"!="") {
			opts_exclusive "vce wboot"    
	}
	if ("`gmm'"!="" & ("`wboot'"!=""|"`wboot1'"!="")) {
		opts_exclusive "gmm wboot"
		local semethod gmm 
	}
	else if ("`gmm'"!="") {
		local semethod "gmm"
	}
	else if (("`wboot'"!=""|"`wboot1'"!="")) {
		if ("`wboot'"!="" & "`wboot1'"!="") {
			display as error "incorrect wildbootstrap specification"
			di as txt "{p 4 4 2}"                           
			di as smcl as err ///
			"You may specify {bf:wboot} or "        			///
			"{bf:wboot()} with arguments but not both."      
			di as smcl as err  "{p_end}"
			exit 198
		}
		local semethod "wildboot"   
		_Parse_Wildboot if `touse', `wboot1' `wboot'
		local seed   "`r(seed)'"
		local reps   = r(reps)
		local wbtype = r(wbtype)
		local cluster "`r(cluster)'"
		if ("`cluster'"!="") {
		    local wncl = r(wncl)
			return scalar wncl = r(wncl)
		}
		return local cluster "`cluster'"
		return local seed "`seed'"
		return local reps = `reps'
		return local wbtype = `wbtype'
	}
	if ("`vce'"!="") {
			_Vce_Parse_Clust, vce(`vce') `gmm'
			return local cluster "`s(cluster)'"
	}
	return local semethod "`semethod'"
end

program _Parse_Wildboot, rclass 
	syntax [anything] [if] [in], [			///
					   WBOOT1				///
					   reps(integer 999) 	///
					   rseed(string) 		///
					   wbtype(string)		///
					   cluster(string)		///
					   ]

	marksample touse 
	
	if ("`wboot1'"=="") {
		return local reps = `reps'
		return local seed "`rseed'"
		if ("`wbtype'"=="") {
		    local wbtypen = 1
		}
		else if ("`wbtype'"=="mammen") {
		    local wbtypen = 1
		}
		else if ("`wbtype'"=="rademacher") {
		    local wbtypen = 2
		}
		else if ("`wbtype'"!="rademacher" & "`wbtype'"!="mammen") {
		    display as error "invalid {bf:wbtype()}"
			di as txt "{p 4 4 2}"                           
			di as smcl as err ///
			"{bf:wbtype()} should be one of {bf:mammen} or " ///
			"{bf:rademacher}."      
			di as smcl as err  "{p_end}"
			exit 198 
		}
		return local wbtype = `wbtypen' 
		if ("`cluster'"!="") {
		    tempvar nclust wncl0
			capture confirm numeric variable `cluster'
			local rc = _rc
			if (`rc') {
				capture destring `rest', generate(`nclust')
				local rc = _rc 
				if (`rc') {
					display in red "option {bf:cluster()} incorrectly specified"
					exit 198
				}
				capture confirm numeric variable `nclust'
				local rc = _rc 
				if (`rc') {
					display in red "option {bf:cluster()} incorrectly specified"
					exit 198
				}
			}
			*quietly egen `wncl0' = group(`cluster') if `touse'
			*summarize `wncl0', meanonly
			*local wncl = r(max)
			*return scalar wncl = `wncl'
		}
		return local cluster "`cluster'"
	}
	else {
	    return local reps = `reps'
		return local seed "`rseed'"
		return local wbtype = 1 
	}
end 


program define _Vce_Parse_Clust, sclass
	syntax [anything], [vce(string) gmm *]
	gettoken key rest : vce, parse(", ")
	local lkey = length(`"`key'"')
	local nvce: list sizeof vce
	local iscluster = 0 
	if (`"`key'"' == bsubstr("cluster",1,max(2,`lkey'))) {
	    local iscluster = 1 
	}
	if (`nvce'>1 & `iscluster'==0 & "`gmm'"=="") {
		display as error "{bf:vce()} option {bf:`key'} not allowed"
		exit 198
	}
	if (`nvce'==1 & "`vce'"!="if" & "`gmm'"=="") {
		display as error "{bf:vce()} option {bf:`key'} not allowed"
		exit 198	    
	}
	if ("`vce'"=="if" & "`gmm'"!="") {
		display as error "{bf:vce()} option {bf:if} not allowed"
		exit 198	    
	}
	if (`nvce'>1) {
		gettoken key rest : vce, parse(", ")
		local lkey = length(`"`key'"')
		local voy = 0 
		if `"`key'"' == bsubstr("cluster",1,max(2,`lkey')) {
			capture confirm numeric variable `rest'
			local rc = _rc
			if (`rc') {
				tempname nclust
				capture destring `rest', generate(`nclust')
				local rc = _rc 
				if (`rc') {
					display in red "option {bf:vce()} incorrectly specified"
					exit 198
				}
				capture confirm numeric variable `nclust'
				local rc = _rc 
				if (`rc') {
					display in red "option {bf:vce()} incorrectly specified"
					exit 198
				}
			}
			local voy = 1
		}
		if ("`key'"=="hac" & "`gmm'"!="") {
			local nvce: list sizeof rest
			local voy = 1
		}
		if ("`key'"=="hac" & "`gmm'"=="") {
			display as error ///
				"{bf:vce()} option {bf:`key'} not allowed with"	///
				" estimator {bf:gmm}"
			exit 198
		}		
		if (`voy'==0) {
				display in red "option {bf:vce()} incorrectly specified"
				exit 198	    
		} 
	}
	if (`iscluster'==1) {
	    local cluster: word 2 of `vce'
	    sreturn local cluster "`cluster'"
	}
end

program define _S_Me_thod, sclass
	if ("`e(method)'"=="drimp") {
		local tmodel "inverse probability tilting"
		local omodel "weighted least squares"
	}
	if ("`e(method)'"=="aipw") {
		local tmodel "inverse probability"
		local omodel "weighted mean"
	}
	if ("`e(method)'"=="dripw") {
		local tmodel "inverse probability"
		local omodel "least squares"
	}
	if ("`e(method)'"=="reg") {
		local tmodel "none"
		local omodel "regression adjustment"
	}
	if ("`e(method)'"=="stdipw") {
		local tmodel "stabilized inverse probability"
		local omodel "weighted mean"
	}
	if ("`e(method)'"=="sipwra") {
		local tmodel "inverse probability"
		local omodel "regression adjustment"
	}

	sreturn local omodel "`omodel'"
	sreturn local tmodel "`tmodel'"
end

program define Display
		syntax [, bmatrix(passthru) vmatrix(passthru) COEFLegend *]
		
        _get_diopts diopts rest, `options'
        local myopts `bmatrix' `vmatrix' 	
		if ("`rest'"!="") {
				display in red "option {bf:`rest'} not allowed"
				exit 198
		}
		
		_S_Me_thod
		local omodel "`s(omodel)'"
		local tmodel "`s(tmodel)'"
		
		if ("`e(method)'"!="all") {
			_coef_table_header, title(Doubly robust difference-in-differences) 
			noi display as text "Outcome model  : {res:`omodel'}"
			noi display as text "Treatment model: {res:`tmodel'}"
			
		}
		else {
			_coef_table_header, ///
				title(Doubly robust difference-in-differences estimator summary) 
		}
		
		if ("`e(policy)'"!="" & "`e(semethod)'"=="wildboot") {
		    if "`e(clustvar)'"!="" {
			    display as text "(Std. err. adjusted for" ///
				as result %9.0gc e(N_clust) ///
				as text " clusters in " as result e(clustvar) as text ")"
			}
			_my_tab_drdid, `diopts'
		}
		if ("`e(semethod)'"!="wildboot") {
		    _coef_table,  `diopts' `myopts' `coeflegend' neq(1)
		}
		
		if ("`e(method)'"=="all") {
			local reg Outcome regression or Regression augmented estimator
			display "{p}Note: This table is provided for comparison" 	///
			" across estimations only. You cannot use it to compare"	///
			" estimates across different estimators{p_end}"
			display "{cmd:dripw} :Doubly Robust IPW"
			display "{cmd:drimp} :Doubly Robust Improved estimator"
			display "{cmd:reg}   :`reg'"
			display "{cmd:ipw}   :Abadie(2005) IPW estimator"
			display "{cmd:stdipw}:Standardized IPW estimator"
			display "{cmd:sipwra}:IPW and Regression adjustment estimator."
		}

end

program _my_tab_drdid, rclass 
	syntax [, level(int `c(level)') noci cformat(string) sformat(string) *]

	_get_diopts diopts rest, `options'

	local cf %9.0g  
	local pf %5.3f
	local sf %7.2f

	if ("`cformat'"!="") {
			local cf `cformat'
	}
	if ("`sformat'"!="") {
			local sf `sformat'
	}

	tempname mytab z t ll ul cimat rtab 
	local policy "`e(policy)'"
	local labn: value label `policy'
	local largo = strlen("`policy'")
	local newvs = 0 
	local novstab  = 0

	if ("`labn'"!="") {
			local uno:  label `labn' 1
			local zero: label `labn' 0 
			local vs "(`uno' vs `zero')"
			local widet = strlen("`vs'")
			local widet0 = `widet'
			local widet = max(`widet', `largo')     
			if (`largo'>`widet0' & `widet'>16) {
					local policy = abbrev("`policy'", 20)
					local widet  = 21 
			}
			if (`widet0'>24) {
					local kn1: list sizeof uno
					local kn0: list sizeof zero
					local kvs0 = max(max(max(`kn1', `kn0'), `largo'), 13)
					local policy = abbrev("`policy'", `kvs0')
					local uno0  = abbrev("`uno'", `kvs0')
					local zero0 = abbrev("`zero'", `kvs0')
					local vs1 "(`uno0' "
					local vs2 "vs"
					local vs3 "`zero0')"
					local newvs = 1 
					local widet = `kvs0' + 2
			}
	}
	else {
			local vs "(1 vs 0)"
			local widet = 13 
			local widet0 = `widet'
			local widet = max(`widet', `largo')
			if (`largo'>`widet0' & `widet'>21) {
					local policy = abbrev("`policy'", 20)
					local widet = 21
			}
	}

	.`mytab' = ._tab.new, col(6) lmargin(0)
	.`mytab'.width    `widet'   |12   12    8     12    12
	.`mytab'.titlefmt  .     .   . %6s    %24s     .
	.`mytab'.pad       .     2   1   0     3     3
	.`mytab'.numfmt    . `cf' `cf' `sf' `cf' `cf'
	
	local stat t 
	/*if "`e(df_r)'" != "" {
			local stat t
			scalar `z' = invttail(e(df_r),(100-`level')/200)
	}
	else {
			local stat z
			scalar `z' = invnorm((100+`level')/200)
	}*/
	local namelist : colname e(b)
	local eqlist : coleq e(b)
	local k : word count `namelist'
	local name : word 1 of `namelist'
	
	matrix `cimat'= e(ciband)
	scalar `ll'   = `cimat'[1,1]
	scalar `ul'   = `cimat'[1,2]
	scalar `t' = `beq'_b[`name']/`beq'_se[`name']
	matrix `rtab' = _b[`name']  \ 	///
	                _se[`name'] \ 	///
					`t'         \ 	///
					.           \	///
					`ll'        \	///
					`ul'        \	///
					.	        \	///
					.	        \	///
					0
	matrix colnames `rtab' = ATET:r1vs0.`policy'
	matrix rownames `rtab' = b se t pvalue ll ul df crit eform
				
	.`mytab'.sep, top
	if `:word count `e(depvar)'' == 1 {
			local depvar "`e(depvar)'"
			local depvar = abbrev("`depvar'", 12)
	}
	.`mytab'.titles "`depvar'"                      /// 1
					" Coefficient"                  /// 2
					"Std. err."						/// 3
					"`stat'"                        /// 4
					"[`level'% conf. interval]" ""  //  5 6

	
	local eq   : word 1 of `eqlist'
	forvalues i = 1/3 {
			if "`eq'" != "_" {
					if "`eq'" != "`eq0'" {
                        .`mytab'.strcolor result  .  .  .  .  .  
                        .`mytab'.strfmt    %-12s  .  .  .  .  .
                        if (`i'==1) {
                                .`mytab'.sep
                                .`mytab'.row      "`eq'" "" "" "" "" ""
                        }
                        .`mytab'.strcolor   text  .  .  .  .  .
                        .`mytab'.strfmt     %12s  .  .  .  .  .

					}
					local beq "[`eq']"
			}
			if (`i'==3) {
			.`mytab'.row    "`vs'"                	///
							`beq'_b[`name']         ///
							`beq'_se[`name']        ///
							`t'                     ///
							`ll' `ul'
			}
			if (`i'==2) {
				.`mytab'.row  "`policy'" "" "" "" "" ""  
			}
	}
	.`mytab'.sep, bottom
	 return matrix table = `rtab'
end

program define easter_egg
		display "{p}This is just for fun. Its my attempt to an Easter Egg within my program. {p_end}" _n /// 
		"{p} Also, if you are reading this, it means you are lucky," ///
		"only 0.1% of people using this program will see this message. {p_end}" _n ///
		"{p} This program was inspired by challenge post by Scott Cunningham. Unforunately, I arrived late due to " ///
		"lack of understanding. I tried to do CSDID (RUN) before learning DRDID (walk). {p_end} " _n  ///
		"{p} Asjad, was the first person who started to properly implement the code in Stata. I got inspired on his work that, " ///
		"Reread the paper, and voala. Everything fit in the form of the first version of this code. {p_end} " _n ///
		"{p} Many thanks go to Pedro, who spend a lot of time explaining details that would otherwise would remain confusing (he is the author fo the original paper after all; {p_end} " _n ///
		"{p} to Miklos, who pushed us to mantain a Github repository for this program. He is taking the side of the user {p_end}" ///
		"{p} To Enrique Pinzon, who helped from the shadows, for a smooth transition from R to Stata (Yay Stata) {p_end}" _n ///
		"{p} and Austin, who was working in parallel when Asjad and I took the lead on this. {p_end}"
end
////////////////////////////////////////////////////////////////////////////////
*** FIRST
**# Abadies
program define drdid_aipw, eclass
	syntax, [						///
			 touse(str) 			///
			 trt(str) 				///		
			 y(str) 				///
			 xvar(str) 				///
			 noisily 				///
			 ivar(str) 				///
			 tmt(str) 				///
			 weight(str) 			///
			 stub(name) 			///	
			 treatvar(string)		///
			 wboot 					///
			 reps(int 999) 			///
			 level(int 95) 			///
			 wbtype(int 1) 			///
			 seed(string)			///
			 cluster(str)			///
			 binit(str)            ///
			 *						///
			 ] 
	** PS
	*set trace on
	_get_diopts diopts other, `options' 
	quietly capture Display, `diopts' `other' 
	
	if _rc==198 {
			Display, `diopts' `other'       
	}
	
	tempvar att psxb __dy__ xb touse2
	tempname psb psV
	
	qui:gen double `att'=.
	if "`ivar'"!="" {
	    *display "Estimating IPW logit"
		qui {		
			`isily' logit `trt' `xvar' if `touse' & `tmt'==0 [iw = `weight']
			predict double `psxb', xb
			
			matrix `psb'=e(b)
			matrix `psV'=e(V)
			** _delta
			bysort `touse' `ivar' (`tmt'):gen double `__dy__'=`y'[2]-`y'[1] if `touse' 
			gen byte `touse2'=`touse'*(`tmt'==0)
  			replace `touse'=0 if `__dy__'==.
 		    tempname b V ciband ncl
			mata:ipw_abadie_panel("`__dy__'","`xvar' ","`xb'","`psb' ","`psV' ","`psxb'","`trt'","`tmt'","`touse2'","`att'","`weight'")
			local ci = `level'/100
			mata:make_tbl("`att'","`cluster' ", "`touse2'", "`b'","`V'","`ciband' ","`ncl'","`wboot' ", `reps', `wbtype', `ci')
			
			matrix colname `b'= ATET:r1vs0.`treatvar'
			matrix colname `V'= ATET:r1vs0.`treatvar' 
			matrix rowname `V'= ATET:r1vs0.`treatvar'
			quietly count if `touse'
			local N = r(N)
			ereturn post `b' `V', buildfvinfo esample(`touse') obs(`N')
			local att1    =`=_b[r1vs0.`treatvar']'
			local attvar1 =`=_se[r1vs0.`treatvar']'^2
			ereturn scalar att1    =`att1'
			ereturn scalar attvar1 =`attvar1'
			ereturn matrix ipwb `psb'
			ereturn matrix ipwV `psV'
			
		}
		*display "DiD with IPW"
		*display "Abadie (2005) inverse probability weighting DiD estimator"
		*ereturn display
	}
	else if "`ivar'"=="" {
	    qui {
			`isily' logit `trt' `xvar' if `touse' [iw = `weight']
			tempvar psxb
			predict double `psxb', xb
			tempname psb psV
			matrix `psb'=e(b)
			matrix `psV'=e(V)
		    tempname b V ciband ncl
 
			mata:ipw_abadie_rc("`y'","`xvar' ","`tmt'","`trt'","`psV'","`psxb'","`weight'","`touse'","`att'")
			** Wbootstrap Multipler
			local ci = `level'/100
			local touse2 `touse'
			mata:make_tbl("`att'","`cluster' ", "`touse2'", "`b'","`V'","`ciband' ","`ncl'","`wboot' ", `reps', `wbtype', `ci')
			
			matrix colname `b' = ATET:r1vs0.`treatvar'
			matrix colname `V' = ATET:r1vs0.`treatvar'
			matrix rowname `V' = ATET:r1vs0.`treatvar'
		}
		
		quietly count if `touse'
		local N = r(N)
		tempvar touse2
		clonevar `touse2'=`touse'
		ereturn post `b' `V', buildfvinfo esample(`touse2') obs(`N')
		local att1    =`=_b[r1vs0.`treatvar']'
		local attvar1 =`=_se[r1vs0.`treatvar']'^2
		ereturn scalar att1    =`att1'
		ereturn scalar attvar1 =`attvar1'
		*ereturn display
		ereturn matrix ipwb `psb'
		ereturn matrix ipwV `psV'
		
	}
	
	
** if STUB is used, then RIF is saved		
	if "`stub'"!="" {
		qui:capture drop `stub'att
		qui:gen double `stub'att=`att'
	}	
	
	if "`cluster'"!="" {
		    ereturn scalar N_clust =`=scalar(`ncl')'
			ereturn local clustvar `cluster'
	}
	if ("`wboot'"=="wboot") {			
		ereturn matrix ciband = `ciband'
		ereturn local semethod "wildboot"
		ereturn local vcetype  "Wboot"
	}
 	ereturn local cmd drdid
	ereturn local method  aipw
	ereturn hidden local method2 aipw
	ereturn local policy "`treatvar'"
	*Display, bmatrix(e(b)) vmatrix(e(V)) `diopts' 
end

** can be more efficient
**#drdid_dripw
program define drdid_dripw, eclass
		syntax, [							///
				touse(str) 					///
				trt(str) 					///
				y(str) 						///
				xvar(str) 					///
				noisily 					///
				ivar(str) 					///
				tmt(str) 					///
				weight(str) 				///
				rc1 						///
				stub(name) 					///
				treatvar(string)			///
				wboot 						///
				level(int 95) 				///
				reps(int 999)				///
				wbtype(int 1) 				///
				seed(string)				///
				cluster(str)				///
				binit(str)            ///
				*							///
				]
	** PS
	tempvar att psxb __dy__ xb touse2
	tempname psb psV regb regV  b V ciband ncl 
	qui:gen double `att'=.
	if "`ivar'"!="" {
	    *display "Estimating IPW logit"
		qui {		
			`isily' logit `trt' `xvar' if `touse' & `tmt'==0 [iw = `weight']
			predict double `psxb', xb			
			matrix `psb'=e(b)
			matrix `psV'=e(V)
			** _delta
			bysort `touse' `ivar' (`tmt'):gen double `__dy__'=`y'[2]-`y'[1] if `touse'
			** Reg for outcome 
 
			`isily' reg `__dy__' `xvar' if `touse' & `trt'==0  & `tmt'==0 [iw = `weight']
			matrix `regb'=e(b)
			matrix `regV'=e(V)
			predict double `xb'
			*capture drop `stub'att
			*gen double `stub'att=. 
			gen byte `touse2'=`touse'*(`tmt'==0)
			replace `touse'=0 if `__dy__'==.
			mata:drdid_panel("`__dy__'","`xvar' ","`xb'","`psb'","`psV'","`psxb'","`trt'","`tmt'","`touse2'","`att'","`weight'")
			**replace `stub'att=. if `tmt'==1
			local ci = `level'/100
			mata:make_tbl("`att'","`cluster' ", "`touse2'", "`b'","`V'","`ciband' ","`ncl'","`wboot' ", `reps', `wbtype', `ci')
			
			matrix colname `b'= ATET:r1vs0.`treatvar'
			matrix colname `V'= ATET:r1vs0.`treatvar'
			matrix rowname `V'= ATET:r1vs0.`treatvar'
			quietly count if `touse'
			local N = r(N)
			ereturn post `b' `V', buildfvinfo esample(`touse') obs(`N')
			local att1    =`=_b[r1vs0.`treatvar']'
			local attvar1 =`=_se[r1vs0.`treatvar']'^2
			ereturn scalar att1    =`att1'
			ereturn scalar attvar1 =`attvar1'
			ereturn matrix ipwb `psb'
			ereturn matrix ipwV `psV'
			ereturn matrix regb `regb'
			ereturn matrix regV `regV'
		}
/*		display "DR DiD with IPW and OLS"
		display "Sant'Anna and Zhao (2020)" _n "{p}Doubly robust DiD estimator based on stabilized inverse probability weighting and ordinary least squares{p_end}"
		ereturn display*/
	}
	else if "`ivar'"=="" {
	    qui {	
			`isily' logit `trt' `xvar' if `touse' [iw = `weight']
			tempvar psxb
			predict double `psxb', xb
			tempname psb psV
			matrix `psb'=e(b)
			matrix `psV'=e(V)
		    tempname b V ciband ncl
			*capture drop `stub'att
			*gen double `stub'att=.
			**ols 
			tempvar y00 y01 y10 y11
			tempname regb00 regb01 regb10 regb11 
			tempname regV00 regV01 regV10 regV11
			`isily' reg `y' `xvar' if `trt'==0 & `tmt' ==0 [iw = `weight']
			predict double `y00'
			matrix `regb00'=e(b)
			matrix `regV00'=e(V)
			`isily' reg `y' `xvar' if `trt'==0 & `tmt' ==1 [iw = `weight']
			predict double `y01'
			matrix `regb01'=e(b)
			matrix `regV01'=e(V)
			`isily' reg `y' `xvar' if `trt'==1 & `tmt' ==0 [iw = `weight']
			predict double `y10'
			matrix `regb10'=e(b)
			matrix `regV10'=e(V)
			`isily' reg `y' `xvar' if `trt'==1 & `tmt' ==1 [iw = `weight']
			predict double `y11'
			matrix `regb11'=e(b)
			matrix `regV11'=e(V)
			if "`rc1'"=="" {
				mata:drdid_rc("`y'","`y00' `y01' `y10' `y11'","`xvar' ","`tmt'","`trt'","`psV'","`psxb'","`weight'","`touse'","`att'")
			}
			else {
			    mata:drdid_rc1("`y'","`y00' `y01' `y10' `y11'","`xvar' ","`tmt'","`trt'","`psV'","`psxb'","`weight'","`touse'","`att'")
				local nle "Not Locally efficient"
			}
			////
			local touse2 `touse'
			local ci = `level'/100
			mata:make_tbl("`att'","`cluster' ", "`touse2'", "`b'","`V'","`ciband' ","`ncl'","`wboot' ", `reps', `wbtype', `ci')			
			matrix colname `b' = ATET:r1vs0.`treatvar'
			matrix colname `V' = ATET:r1vs0.`treatvar'
			matrix rowname `V' = ATET:r1vs0.`treatvar'			
		}
		*display "DR DiD with IPW and OLS for Repeated Crossection: `nle'"
		*display "Sant'Anna and Zhao (2020)" _n "{p}Doubly robust DiD estimator based on stabilized inverse probability weighting and ordinary least squares{p_end}"
		quietly count if `touse'
		local N = r(N)
		tempvar touse2	
		clonevar `touse2'=`touse'
		ereturn post `b' `V', buildfvinfo esample(`touse2') obs(`N')
		*ereturn post `b' `V'
		*ereturn display
		local att1    =`=_b[r1vs0.`treatvar']'
		local attvar1 =`=_se[r1vs0.`treatvar']'^2
		ereturn scalar att1    =`att1'
		ereturn scalar attvar1 =`attvar1'
 		ereturn matrix ipwb `psb'
		ereturn matrix ipwV `psV'
		
		ereturn matrix regb00 `regb00'
		ereturn matrix regV00 `regV00'
		ereturn matrix regb01 `regb01'
		ereturn matrix regV01 `regV01'
		ereturn matrix regb10 `regb10'
		ereturn matrix regV10 `regV10'
		ereturn matrix regb11 `regb11'
		ereturn matrix regV11 `regV11'
 
	}
	if "`stub'"!="" {
		qui:capture drop `stub'att
		qui:gen double `stub'att=`att'
	}
	
	if "`cluster'"!="" {
		    ereturn scalar N_clust =`=scalar(`ncl')'
			ereturn local clustvar `cluster'
	}
 
	if ("`wboot'"=="wboot") {
		ereturn matrix ciband = `ciband'
		ereturn local semethod "wildboot"
	}
	
	ereturn local cmd drdid
	ereturn local method  dripw
	ereturn hidden local method2 dripw
	ereturn local policy "`treatvar'"
	*Display, bmatrix(e(b)) vmatrix(e(V)) `diopts' 
end

**#drdid_reg
program define drdid_reg, eclass
	syntax, [						///
			touse(str) 				///
			trt(str) 				///
			y(str) 					///
			xvar(str) 				///
			noisily 				///
			ivar(str) 				///
			tmt(str) 				///
			weight(str) 			///
			stub(name) 				///
			treatvar(string) 		///
			wboot 					///
			level(int 95) 			///
			reps(int 999) 			///
			wbtype(int 1) 			///
			seed(string)			///
			 cluster(str)			///
			 binit(str)            ///
			*						///
			] 
** Simple application. But right now without RIF
	tempvar att
	qui:gen double `att'=.
	if "`ivar'"!="" {
		qui {
			tempvar __dy__ xb touse2
			tempname regb regV b V ciband ncl
			bysort `touse' `ivar' (`tmt'):gen double	///
				`__dy__'=`y'[2]-`y'[1] if `touse' 
				`isily' reg `__dy__' `xvar' if `touse' & `trt'==0	///
				& `tmt'==0  [iw = `weight']
			predict double `xb'
			//////////////////////		
			matrix `regb'=e(b)
			matrix `regV'=e(V)
			gen byte `touse2'=`touse'*(`tmt'==0)
			replace `touse'=0 if `__dy__'==.
			mata:reg_panel("`__dy__'", "`xvar' ", "`xb' " , "`trt'",	///
				"`tmt'" , "`touse2'","`att'","`weight'") 
			local ci = `level'/100
			mata:make_tbl("`att'","`cluster' ", "`touse2'", "`b'","`V'","`ciband' ","`ncl'","`wboot' ", `reps', `wbtype', `ci')			
			matrix colname `b' = ATET:r1vs0.`treatvar'
			matrix colname `V' = ATET:r1vs0.`treatvar'
			matrix rowname `V' = ATET:r1vs0.`treatvar'
			quietly count if `touse'
			local N = r(N)
			ereturn post `b' `V', buildfvinfo esample(`touse') obs(`N')
			local att1    =`=_b[r1vs0.`treatvar']'
			local attvar1 =`=_se[r1vs0.`treatvar']'^2
			ereturn scalar att1    =`att1'
			ereturn scalar attvar1 =`attvar1'
			ereturn matrix regb `regb'
			ereturn matrix regV `regV' 
		}
		*display "DiD with OR" _n "Outcome regression DiD estimator based on ordinary least squares"
		*ereturn display
	}	
	else if "`ivar'"=="" {
		qui {
		    tempvar y00 y01
			tempname regb00 regb01 regV00 regV01
		    `isily' reg `y' `xvar' if `trt'==0 & `tmt' ==0 [iw = `weight']
			predict double `y00'
			matrix `regb00' = e(b)
			matrix `regV00' = e(V)
			`isily' reg `y' `xvar' if `trt'==0 & `tmt' ==1 [iw = `weight']
			predict double `y01'
			matrix `regb01' = e(b)
			matrix `regV01' = e(V)
			tempname b V ciband ncl
			*capture drop `stub'att
			*gen double `stub'att=.
			mata:reg_rc("`y'","`y00' `y01'","`xvar' ",	///
				"`tmt'","`trt'","`weight'","`touse'","`att'")
			local ci = `level'/100
			local touse2 `touse'
			mata:make_tbl("`att'","`cluster' ", "`touse2'", "`b'","`V'","`ciband' ","`ncl'","`wboot' ", `reps', `wbtype', `ci')			
			matrix colname `b' = ATET:r1vs0.`treatvar'
			matrix colname `V' = ATET:r1vs0.`treatvar'
			matrix rowname `V' = ATET:r1vs0.`treatvar'
			quietly count if `touse'
			local N = r(N)
			tempvar touse2
			clonevar `touse2'=`touse'
			ereturn post `b' `V', buildfvinfo esample(`touse2') obs(`N')
			ereturn matrix regb00 `regb00'
			ereturn matrix regV00 `regV00' 
			ereturn matrix regb01 `regb01'
			ereturn matrix regV01 `regV01'
			local att1    =`=_b[r1vs0.`treatvar']'
			local attvar1 =`=_se[r1vs0.`treatvar']'^2
			ereturn scalar att1    =`att1'
			ereturn scalar attvar1 =`attvar1'
		}
		*display "DiD with OR for RC" _n "Outcome regression DiD estimator based on ordinary least squares"
		*ereturn display
	}
	
	if "`stub'"!="" {
		qui:capture drop `stub'att
		qui:gen double `stub'att=`att'
	}
	if "`cluster'"!="" {
		    ereturn scalar N_clust =`=scalar(`ncl')'
			ereturn local clustvar `cluster'
		}
		
	if ("`wboot'"=="wboot") {			
		ereturn matrix ciband = `ciband'
		ereturn local semethod "wildboot"
	}
	ereturn local cmd drdid
	ereturn local method  reg
	ereturn hidden local method2 reg
	ereturn local policy "`treatvar'"
	*Display, bmatrix(e(b)) vmatrix(e(V)) `diopts' 
	
end

**#drdid_sipw
program define drdid_stdipw, eclass
	syntax, [						///
			touse(str) 				///
			trt(str) 				///
			y(str) 					///
			xvar(str) 				///
			noisily 				///
			ivar(str) 				///
			tmt(str) 				///
			weight(str) 			///
			stub(name) 				///
			treatvar(string) 		///
			wboot 					///
			level(int 95) 			///
			reps(int 999) 			///
			wbtype(int 1) 			///
			seed(string)			///
			 cluster(str)			///
			 binit(str)            ///
			*						///
			] 
			
	tempvar att
	qui:gen double `att'=.
** Simple application. But right now without RIF
	if "`ivar'"!="" {
	    *display "Estimating IPW logit"
		qui {		
			`isily' logit `trt' `xvar' if `touse' & `tmt'==0 [iw = `weight']
			tempvar psxb __dy__ xb touse2
			tempname psb psV b V ciband ncl
			predict double `psxb', xb		
			matrix `psb'=e(b)
			matrix `psV'=e(V)
			** _delta
			bysort `touse' `ivar' (`tmt'):gen double `__dy__'=`y'[2]-`y'[1]	///
				if `touse'
			** Reg for outcome 
		}
		*display "Estimating Counterfactual Outcome"	
		qui {
			*`isily' reg `__dy__' `xvar' if `trt'==0 
			*gen double `stub'att=.		
			gen byte `touse2'=`touse'*(`tmt'==0)
			replace `touse'=0 if `__dy__'==.
			mata:std_ipw_panel("`__dy__'","`xvar' ","`xb'",		///
				"`psb'","`psV'","`psxb'","`trt'","`tmt'","`touse2'",	///
				"`att'","`weight'")
			local ci = `level'/100
			mata:make_tbl("`att'","`cluster' ", "`touse2'", "`b'","`V'","`ciband' ","`ncl'","`wboot' ", `reps', `wbtype', `ci')			
			matrix colname `b'= ATET:r1vs0.`treatvar'
			matrix colname `V'= ATET:r1vs0.`treatvar'
			matrix rowname `V'= ATET:r1vs0.`treatvar'

		}
		*display "DiD with stabilized IPW" _n "{p}Abadie (2005) inverse probability weighting DiD estimator with stabilized weights{p_end}" 
		quietly count if `touse'
		local N = r(N)
		ereturn post `b' `V', buildfvinfo esample(`touse') obs(`N')
		*ereturn display
		local att1    =`=_b[r1vs0.`treatvar']'
		local attvar1 =`=_se[r1vs0.`treatvar']'^2
		ereturn scalar att1    =`att1'
		ereturn scalar attvar1 =`attvar1'
		ereturn matrix ipwb `psb'
		ereturn matrix ipwV `psV'
	}	
	else if "`ivar'"=="" {
	    qui {
		    `isily' logit `trt' `xvar' if `touse' [iw = `weight']
			tempvar psxb
			predict double `psxb', xb
			tempname psb psV
			matrix `psb'=e(b)
			matrix `psV'=e(V)
			tempname b V ciband ncl
			*capture drop `stub'att
			*gen `stub'att=.
			mata:std_ipw_rc("`y'","`xvar' ","`tmt'","`trt'","`psV'","`psxb'","`weight'","`touse'","`att'")
			local ci = `level'/100
			local touse2 `touse'
			mata:make_tbl("`att'","`cluster' ", "`touse2'", "`b'","`V'","`ciband' ","`ncl'","`wboot' ", `reps', `wbtype', `ci')			
			matrix colname `b' = ATET:r1vs0.`treatvar'
			matrix colname `V' = ATET:r1vs0.`treatvar'
			matrix rowname `V' = ATET:r1vs0.`treatvar'
		}
		
		quietly count if `touse'
		local N = r(N)
		tempvar touse2
		clonevar `touse2'=`touse'
		ereturn post `b' `V', buildfvinfo esample(`touse2') obs(`N') 
		local att1    =`=_b[r1vs0.`treatvar']'
		local attvar1 =`=_se[r1vs0.`treatvar']'^2
		ereturn scalar att1    =`att1'
		ereturn scalar attvar1 =`attvar1'
		ereturn matrix ipwb `psb'
		ereturn matrix ipwV `psV'
		if "`cluster'"!="" {
		    ereturn scalar N_clust =`=scalar(`ncl')'
			ereturn local clustvar `cluster'
		}
	}
	
	if "`stub'"!="" {
		qui:capture drop `stub'att
		qui:gen double `stub'att=`att'
	}
	if "`cluster'"!="" {
		    ereturn scalar N_clust =`=scalar(`ncl')'
			ereturn local clustvar `cluster'
		}
	if ("`wboot'"=="wboot") {			
		ereturn matrix ciband = `ciband'
		ereturn local semethod "wildboot"
	}
	ereturn local cmd drdid
	ereturn local method stdipw 
	ereturn hidden local method2 stdipw 
	ereturn local policy "`treatvar'"
	*Display, bmatrix(e(b)) vmatrix(e(V)) `diopts' 
	
end

// only one without Mata writting. Consider working on it
**#drdid_sipwra
program define drdid_sipwra, eclass
	syntax, [						///
			touse(str) 				///
			trt(str) 				///
			y(str) 					///
			xvar(str) 				///
			noisily 				///
			ivar(str) 				///
			tmt(str) 				///
			weight(str) 			///
			stub(name) 				///
			treatvar(string)		///
			reps(int 999) 			/// Notused here
			wbtype(int 1) 			/// not used here
			seed(string)			///
			level(int 95) 			/// not used here
			cluster(str)			///
			*						///
			] 
	tempvar att
	qui:gen double `att'=.
** Simple application. But right now without RIF
	if "`cluster'"!="" {
		local clopt vce(cluster `cluster')
	}
   if "`ivar'"=="" {
       display "Estimator not implemented for RC"
	   exit 198 
   }
   else {
	qui {
		tempvar __dy__ sy
		bysort `touse' `ivar' (`tmt'):gen double `__dy__'=`y'[2]-`y'[1] if `touse'
		sum `__dy__' if  `touse'
		local scl = r(mean)
		gen double `sy'= `__dy__'/`scl'		
		qui:teffects ipwra (`sy' `xvar') (`trt' `xvar', logit)	///
			if `touse' & `tmt'==0 [iw = `weight'] , atet `clopt' iter(5)
		tempname b V ciband ncl aux
		matrix `aux'=e(b)*`scl'
		matrix `b'=`aux'[1,1]
		matrix `aux'=e(V)*`scl'^2
		matrix `V'=`aux'[1,1]
		matrix colname `b' = ATET:r1vs0.`treatvar'
		matrix colname `V' = ATET:r1vs0.`treatvar'
		matrix rowname `V' = ATET:r1vs0.`treatvar'
		quietly count if `touse'
		local N = r(N)
		ereturn post `b' `V', buildfvinfo esample(`touse') obs(`N')
	}
		*ereturn display	
   }

	if "`stub'"!="" {
		qui:capture drop `stub'att
		qui:gen double `stub'att=`att'
	}   
	ereturn local cmd drdid
	ereturn local method sipwra
	ereturn hidden local method2 sipwra
	ereturn local policy "`treatvar'"
	*Display, bmatrix(e(b)) vmatrix(e(V)) `diopts' 
end
 
**#drdid_dript
program define drdid_imp, eclass  
	syntax, [					///
			touse(str)			///
			trt(str)			///
			y(str) 				///
			xvar(str) 			///
			noisily 			///
			ivar(str) 			///
			tmt(str) 			///
			weight(str) 		///
			rc1 stub(name)		///
			treatvar(string)	///
            wboot 				///
			level(int 95) 		///
			reps(int 999) 		///
			wbtype(int 1) 		///
			seed(string)			///
			cluster(str)			///
			binit(str)            ///
		*						///
		] 
 
	_get_diopts diopts other, `options' 
	quietly capture Display, `diopts' `other' 
	
	if _rc==198 {
			Display, `diopts' `other'       
	}
	
	tempvar att  psxb __dy__ w0 xb
	tempname iptb iptV regb regV b V ciband ncl

	qui:gen double `att'=.
	
	if "`ivar'"!="" {
		qui {
			
			`isily'  mlexp (`trt'*{xb:`xvar' _cons}-(`trt'==0)*exp({xb:}))  ///
					if `touse' & `tmt'==0 [iw = `weight'],  from(`binit') ///
					 derivative(/xb=`trt'-(`trt'==0)*exp({xb:}))
			*vce(robust)
			matrix `iptb'=e(b)
			matrix `iptV'=e(V)
			predict double `psxb',xb
			
			** Determine dy and dyhat
			bysort `touse' `ivar' (`tmt'):gen double `__dy__'=`y'[2]-`y'[1] if `touse'

			** determine weights
			gen double `w0' = ((logistic(`psxb')*(1-`trt')))/(1-logistic(`psxb'))*`weight'
			sum `w0' if `touse' , meanonly
			replace `w0'=`w0'/r(mean)
			
			** estimating dy_hat for a counterfactual
			`isily' reg `__dy__' `xvar' [w=`w0'] if `trt'==0 & `tmt'==0,
			matrix `regb' =e(b)
			matrix `regV' =e(V)
			qui:predict double `xb'
			tempvar touse2
			gen byte `touse2'=`touse'*(`tmt'==0)
			replace `touse'=0 if `__dy__'==.
			mata:drdid_imp_panel("`__dy__'","`xvar' ","`xb'",	///
				"`psb'","`psV'","`psxb'","`trt'","`tmt'","`touse2'",	///
				"`att'","`weight'")	
********************************************************************************				
********************************************************************************
			
			local ci = `level'/100
			mata:make_tbl("`att'","`cluster' ", "`touse2'", "`b'","`V'","`ciband' ","`ncl'","`wboot' ", `reps', `wbtype', `ci')
			
 			matrix colname `b'= ATET:r1vs0.`treatvar'
			matrix colname `V'= ATET:r1vs0.`treatvar'
			matrix rowname `V'= ATET:r1vs0.`treatvar'
			quietly count if `touse'
			local N = r(N)
			ereturn post `b' `V', buildfvinfo esample(`touse') obs(`N')
			local att1    =`=_b[r1vs0.`treatvar']'
			local attvar1 =`=_se[r1vs0.`treatvar']'^2
			ereturn scalar att1    =`att1'
			ereturn scalar attvar1 =`attvar1'
			ereturn matrix iptb `iptb'
			ereturn matrix iptV `iptV'
			ereturn matrix regb `regb'
			ereturn matrix regV `regV'
			ereturn local cmd drdid
			ereturn local method drimp
 
			if ("`wboot'"=="wboot") {			
				ereturn matrix ciband = `ciband'
				ereturn local semethod "wildboot"
			}
		}
		*display "DR DiD with IPT and WLS" _n "{p}Sant'Anna and Zhao (2020) Improved doubly robust DiD estimator based on inverse probability of tilting and weighted least squares{p_end}"
		*ereturn display
	}
	else {
	**# for Crossection estimator				

		qui {

			/*`isily' gmm ((`trt'==1)-(`trt'==0)*exp({b:`xvar' _cons})) if `touse'  [iw = `weight'], ///
			instrument(`xvar' ) derivative(/b=-(`trt'==0)*exp({b:})) ///
			onestep winit(identity) */
			
			`isily'  mlexp (`trt'*{xb:`xvar' _cons}-(`trt'==0)*exp({xb:}))  ///
					if `touse' [iw = `weight'], vce(robust) from(`binit') ///
					 derivative(/xb=`trt'-(`trt'==0)*exp({xb:}))

			//& `tmt'==0 
			tempname iptb iptV regb00 regV00 regb01 regV01 regb10	///
					 regV10 regb11 regV11
			tempvar psxb w1 w0 y01 y00 y10 y11
			
			matrix `iptb'=e(b)
			matrix `iptV'=e(V)
			predict double `psxb', xb
			** outcomes
			gen double `w0' = (1-`trt')*logistic(`psxb')/(1-logistic(`psxb'))
			`isily' reg `y' `xvar' [w=`w0'] if `trt'==0 & `tmt'==0,
			predict double `y00'
			matrix `regb00' =e(b)
			matrix `regV00' =e(V)
			`isily' reg `y' `xvar' [w=`w0'] if `trt'==0 & `tmt'==1,
			predict double `y01'
			matrix `regb01' =e(b)
			matrix `regV01' =e(V)
			`isily' reg `y' `xvar'  		   if `trt'==1 & `tmt'==0,
			predict double `y10'
			matrix `regb10' =e(b)
			matrix `regV10' =e(V)
			`isily' reg `y' `xvar'  		   if `trt'==1 & `tmt'==1,
			predict double `y11'
			matrix `regb11' =e(b)
			matrix `regV11' =e(V)
			tempname b V ciband ncl
			if "`rc1'"=="" {
				mata:drdid_imp_rc("`y'","`y00' `y01' `y10' `y11'",	///
					"`xvar' ","`tmt'","`trt'","`iptV'","`psxb'",	///
					"`weight'","`touse'","`att'")
			}
			else {
			    mata:drdid_imp_rc1("`y'","`y00' `y01' `y10' `y11'",	///
					"`xvar' ","`tmt'","`trt'","`iptV'","`psxb'",	///
					"`weight'","`touse'","`att'")
				local nle "Not Locally efficient"
			}
			local ci = `level'/100
			local touse2 `touse'
			mata:make_tbl("`att'","`cluster' ", "`touse2'", "`b'","`V'","`ciband' ","`ncl'","`wboot' ", `reps', `wbtype', `ci')
			  			
			matrix colname `b'=ATET:r1vs0.`treatvar'
			matrix colname `V'=ATET:r1vs0.`treatvar'
			matrix rowname `V'=ATET:r1vs0.`treatvar'
			quietly count if `touse'
			local N = r(N)
			tempvar touse2	
			clonevar `touse2'=`touse'
			ereturn post `b' `V', buildfvinfo esample(`touse2') obs(`N')
			ereturn local cmd drdid
			ereturn local method drimp
			
			if ("`wboot'"=="wboot") {			
				ereturn matrix ciband = `ciband'
				ereturn local semethod "wildboot"
			}
		}
		
		*display "{p}DR DiD with IPT and WLS for OLS `nle'{p_end}" _n "{p}Sant'Anna and Zhao (2020) Improved doubly robust DiD estimator based on inverse probability of tilting and weighted least squares{p_end}"
		*ereturn display
		local att1    =`=_b[r1vs0.`treatvar']'
		local attvar1 =`=_se[r1vs0.`treatvar']'^2
		ereturn scalar att1    =`att1'
		ereturn scalar attvar1 =`attvar1'
		ereturn matrix iptb `iptb'
		ereturn matrix iptV `iptV'
		ereturn matrix regb00 `regb00'
		ereturn matrix regV00 `regV00'
		ereturn matrix regb01 `regb01'
		ereturn matrix regV01 `regV01'
		ereturn matrix regb10 `regb10'
		ereturn matrix regV10 `regV10'
		ereturn matrix regb11 `regb11'
		ereturn matrix regV11 `regV11'
	}
	if "`stub'"!="" {
		qui:capture drop `stub'att
		qui:gen double `stub'att=`att'
	}
	if "`cluster'"!="" {
		    ereturn scalar N_clust =`=scalar(`ncl')'
			ereturn local clustvar `cluster'
	}
	ereturn local policy "`treatvar'"
	*Display, bmatrix(e(b)) vmatrix(e(V)) `diopts' 
end

program define _S__est, sclass
	syntax [anything], [drimp dripw reg stdipw ipw ipwra all * ]
	
	if ("`drimp'"!="") {
		local drimp imp 
	}
	if ("`ipw'"!="") {
		local ipw "aipw"
	}
	if ("`ipwra'"!="") {
		local ipwra "sipwra"
	}
	sreturn local estimator ///
		"`drimp'`dripw'`reg'`stdipw'`ipw'`ipwra'`all'"
end 

mata 
	////////////////////////////////////////////////////////////////////////////////////////////////
// dript
	void drdid_imp_panel(string scalar dy_, xvar_, xb_ , psb_,psV_,psxb_,trt_,tmt_,touse,rif,ww) {
	    real matrix dy, xvar, xb, psb, psv, psxb, trt, tmt
		// This code is based on Asjad Replication
		// Gather all data
		real scalar nn
		dy  =st_data(.,dy_  ,touse)
		nn=rows(dy)
		//st_view(xvar=.,xvar_,touse)
		xb=st_data(.,xb_,touse)
		psxb=st_data(.,psxb_,touse)
		real matrix psc
		psc=logistic(psxb)
		trt =st_data(.,trt_ ,touse)
		real matrix w
		w=st_data(.,ww,touse)
		real matrix w_1 , w_0, att, att_inf_func
		w_1 = w :* trt
		w_0 = w :* psc :* (1:-trt):/(1:-psc)
		w_1 = w_1:/mean(w_1)
		w_0 = w_0:/mean(w_0)
		att=(dy:-xb):*(w_1:-w_0)
		att_inf_func = mean(att) :+ att :- w_1:*mean(att)
		st_store(.,rif,touse, att_inf_func)	
		// Variance should be divided by n not n-1
		
		
	}
 
	// standard IPW
  	void std_ipw_panel(string scalar dy_, xvar_, xb_ , psb_,psV_,psxb_,trt_,tmt_,touse,rif,ww) {
	    real matrix dy, xvar, xb, psb, psv, psxb, trt, tmt
 		// Gather all data
		real scalar nn
		dy  =st_data(.,dy_  ,touse)
		nn=rows(dy)
		// verify xvar
		if (xvar_==" ") {
			xvar=J(rows(dy),1,1)	
		}
		else xvar=st_data(.,xvar_,touse),J(rows(dy),1,1)	
		
		//xb=st_data(.,xb_,touse)
		psxb=st_data(.,psxb_,touse)
		real matrix psc
		psc=logistic(psxb)
		trt =st_data(.,trt_ ,touse)
		tmt =st_data(.,tmt_ ,touse)
		// and matrices
		//psb =st_matrix(psb_ )
		psv =st_matrix(psV_ )
		// for now assume weights = 1
		real matrix w
		w=st_data(.,ww,touse)
		
		real matrix w_1, w_0, att_cont, att_treat,
					eta_treat, eta_cont, 
					lin_ps,  att_inf_func
			
		w_1= w :* trt
		w_0= w :* psc :* (1 :- trt):/(1 :- psc)
		att_treat = w_1:* dy
		att_cont  = w_0:* dy
		eta_treat = mean(att_treat)/mean(w_1)
		eta_cont  = mean(att_cont)/mean(w_0)
		ipw_att   = eta_treat :- eta_cont
		
		inf_treat = (att_treat :- w_1 :* eta_treat)/mean(w_1)
		inf_cont_1 = (att_cont :- (w_0 :* eta_cont))
		lin_ps = (w:* (trt :- psc) :* xvar)*(psv * nn)
		//M2 =
		inf_cont_2 = lin_ps * mean(w_0 :* (dy :- eta_cont) :* xvar)'
		inf_control = (inf_cont_1 :+ inf_cont_2)/mean(w_0)
		att_inf_func = mean(ipw_att):+inf_treat :- inf_control
 
		st_store(.,rif,touse, att_inf_func)	

	}
  
	void reg_panel(string scalar dy_, xvar_, xb_ , trt_,tmt_,touse,rif,ww) {
	    real matrix dy, xvar, xb, trt, tmt
		// This code is based on Asjad Replication
		// Gather all data
		real scalar nn
		dy  =st_data(.,dy_  ,touse)
		nn=rows(dy)
		
		// verify xvar
		if (xvar_==" ") {
			xvar=J(rows(dy),1,1)	
		}
		else xvar=st_data(.,xvar_,touse),J(rows(dy),1,1)	
		
		xb=st_data(.,xb_,touse)
		// psxb=st_data(.,psxb_,touse)
		// real matrix psc
		// psc=logistic(psxb)
		trt =st_data(.,trt_ ,touse)
		tmt =st_data(.,tmt_ ,touse)
		// and matrices
		// psb =st_matrix(psb_ )
		// psv =st_matrix(psV_ )
		// for now assume weights = 1
		real matrix w
		w=st_data(.,ww,touse)
		real matrix w_1, w_0, att_cont, att_treat,
					eta_treat, eta_cont, wols_x, wols_eX,
					lin_ols, att_inf_func
		
		w_1 = w :* trt
		w_0 = w :* trt
		att_treat = w_1:* dy
		att_cont  = w_0:* xb
		eta_treat = mean(att_treat):/mean(w_1)
		eta_cont  = mean(att_cont)  :/mean(w_0)
		reg_att = eta_treat :- eta_cont
		wols    = w :* (1 :- trt)
		wols_x  = wols :* xvar
		wols_eX = wols :* (dy:-xb) :* xvar
		
		XpX_inv = invsym(quadcross(wols_x, xvar))*nn
		lin_ols = wols_eX * XpX_inv
		inf_treat    = (att_treat :- w_1 * eta_treat):/mean(w_1)
		inf_cont_1   = (att_cont :- w_0 * eta_cont)
		inf_cont_2   = lin_ols * mean(w_0 :* xvar )'
		inf_control  = (inf_cont_1 :+ inf_cont_2):/mean(w_0)
		att_inf_func = mean(reg_att):+(inf_treat :- inf_control)
			
		st_store(.,rif,touse, att_inf_func)	

				
	}
	
 
/////////////////////////////////////////////////////////////////////////////////////////////////	
// TIPW drdid_tipw
 	void ipw_abadie_panel(string scalar dy_, xvar_, xb_ , psb_,psV_,psxb_,trt_,tmt_,touse,rif,ww) {
	    real matrix dy, xvar, xb, psb, psv, psxb, trt, tmt
 		// Gather all data
		real scalar nn
		dy  =st_data(.,dy_  ,touse)
		nn=rows(dy)
		// verify xvar
		if (xvar_==" ") {
			xvar=J(rows(dy),1,1)	
		}
		else xvar=st_data(.,xvar_,touse),J(rows(dy),1,1)	
		
		//xb=st_data(.,xb_,touse)
		psxb=st_data(.,psxb_,touse)
		real matrix psc
		psc=logistic(psxb)
		trt =st_data(.,trt_ ,touse)
		tmt =st_data(.,tmt_ ,touse)
		// and matrices
		psb =st_matrix(psb_ )
		psv =st_matrix(psV_ )
		// for now assume weights = 1
		real matrix w
		w=st_data(.,ww,touse)
		real matrix w_1, w_0, att_cont, att_treat,
					eta_treat, eta_cont, ipw_att,
					lin_ps, att_lin1, mom_logit, att_lin2, att_inf_func
					
		w_1= w :* trt
		w_0= w :* psc :* (1 :- trt):/(1 :- psc)
		att_treat = w_1:* dy
		att_cont  = w_0:* dy
		eta_treat = mean(att_treat)/mean(w_1)
		eta_cont  = mean(att_cont)/mean(w_1)
		ipw_att   = eta_treat - eta_cont
		lin_ps = (w:* (trt :- psc) :* xvar)*(psv * nn)
		att_lin1  = att_treat :- att_cont
		mom_logit = mean(att_cont  :* xvar)
		att_lin2  = lin_ps * mom_logit'
		att_inf_func = mean(ipw_att):+(att_lin1 :- att_lin2 :- w_1 :* ipw_att)/mean(w_1)
 		st_store(.,rif,touse, att_inf_func)	
		// Variance should be divided by n not n-1

		
	}
	
	
/// drdid_ipw	
 	void drdid_panel(string scalar dy_, xvar_, xb_ , psb_,psV_,psxb_,trt_,tmt_,touse,rif,ww) {
	    real matrix dy, xvar, xb, psb, psv, psxb, trt, tmt
		// This code is based on Asjad Replication
		// Gather all data
		real scalar nn
		dy  =st_data(.,dy_  ,touse)
		nn=rows(dy)
		// verify xvar
		if (xvar_==" ") {
			xvar=J(rows(dy),1,1)	
		}
		else xvar=st_data(.,xvar_,touse),J(rows(dy),1,1)			
		//xvar=st_data(.,xvar_,touse),J(rows(dy),1,1)
		xb=st_data(.,xb_,touse)
		psxb=st_data(.,psxb_,touse)
		real matrix psc
		psc=logistic(psxb)
		trt =st_data(.,trt_ ,touse)
		tmt =st_data(.,tmt_ ,touse)
		// and matrices
		psb =st_matrix(psb_ )
		psv =st_matrix(psV_ )
		// for now assume weights = 1
		real matrix w
		w=st_data(.,ww,touse)
		
		// TRAD DRDID
		real matrix w_1, w_0 
		w_1 = (w:*trt)
		w_1 = w_1 /mean(w_1)
		w_0 = (w:*psc:*(-trt:+1):/(-psc:+1))
		w_0 = w_0 /mean(w_0 ) 
		// ATT
		real matrix dy_xb, att, w_ols, wols_eX,
					XpX_inv, lin_wols, lin_ps, 
					n1, n0, nest, a, att_inf_func
					
		dy_xb=dy:-xb
		att = mean((w_1:-w_0):*(dy_xb))
		
		// influence functions OLS
		w_ols 	 =    w :* (1 :- trt)
		wols_eX  = w_ols:* (dy_xb):* xvar
		XpX_inv  = invsym(quadcross(xvar,w_ols,xvar))*nn 
		lin_wols =  wols_eX * XpX_inv   
		// IF for logit
		lin_ps 	   = (w :* (trt:-psc) :* xvar) * (psv * nn)
		// Components for RIF
		n1   = w_1:*((dy_xb):-mean(dy_xb,w_1))
		n0   = w_0:*((dy_xb):-mean(dy_xb,w_0))
		
		a    = ((1:-trt):/(1:-psc):^2)/ mean(psc:*(1:-trt):/(1:-psc))
		// This only works because w_1 and w_0 are mutually exclusive
		nest = lin_wols * (mean(xvar,w_1):-mean(xvar,w_0))' :+
		       lin_ps   * mean( a :* (dy_xb :- mean(dy_xb,w_0)) :* exp(psxb):/(1:+exp(psxb)):^2:*xvar)'			   
		// RIF att_inf_func = inf_treat' :- inf_control
		att_inf_func = att:+n1:-n0:-nest
		st_store(.,rif,touse, att_inf_func)	
		// Variance should be divided by n not n-1
	
	}
	//// standard IPW

//# RC
/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////	
	 void std_ipw_rc(string scalar y_, xvar_ , tmt_, trt_, psv_, pxb_, wgt_ ,  touse, rif){
    // main Loading variables
    real matrix y, 	xvar, tmt, trt, psv, psc, wgt
				
	y    = st_data(.,y_, touse)
		// verify xvar
		if (xvar_==" ") {
			xvar=J(rows(y),1,1)	
		}
		else xvar=st_data(.,xvar_,touse),J(rows(y),1,1)		
	//xvar = st_data(.,xvar_, touse), J(rows(y),1,1)
	tmt  = st_data(.,tmt_, touse)
	trt  = st_data(.,trt_, touse)
	psc  = logistic(st_data(.,pxb_, touse))
	wgt  = st_data(.,wgt_, touse)
	psv  = st_matrix(psv_)
	real scalar nn
 	nn = rows(y)
	
	real matrix w10, w11, w00, w01
	
    w10 = wgt :* trt :* (1 :- tmt)
    w11 = wgt :* trt :* tmt
    w00 = wgt :* psc :* (1 :- trt) :* (1 :- tmt):/(1 :- psc)
    w01 = wgt :* psc :* (1 :- trt) :* tmt:/(1 :- psc)
	
	w00 = w00:/mean(w00 )
	w01 = w01:/mean(w01 )
	w10 = w10:/mean(w10 )
	w11 = w11:/mean(w11 )
	
	real matrix att_treat_pre, att_treat_post, att_cont_pre, att_cont_post,
				eta_treat_pre, eta_treat_post, eta_cont_pre, eta_cont_post
	att_treat_pre  	= w10 :* y
    att_treat_post 	= w11 :* y
    att_cont_pre   	= w00 :* y
    att_cont_post  	= w01 :* y
    eta_treat_pre  	= mean(att_treat_pre)
    eta_treat_post 	= mean(att_treat_post)
    eta_cont_pre   	= mean(att_cont_pre)
    eta_cont_post  	= mean(att_cont_post)
	
	real matrix ipw_att, lin_rep_ps
    ipw_att 		= (eta_treat_post :- eta_treat_pre) :- (eta_cont_post :- eta_cont_pre)
    //score_ps 		= wgt :* (trt :- psc) :* xvar
    //Hessian_ps 		= psv :* nn
    lin_rep_ps 		= (wgt :* (trt :- psc) :* xvar) * (psv :* nn)
    
	real matrix inf_treat, inf_cont, M2_pre, M2_post, inf_cont_ps, att_inf_func
	
	inf_treat 		= (att_treat_post :- w11 :* eta_treat_post) :- (att_treat_pre :- w10 :* eta_treat_pre)
    inf_cont 		= (att_cont_post :- w01 :* eta_cont_post) :- (att_cont_pre :- w00 :* eta_cont_pre)
    M2_pre 			= mean(w00 :* (y :- eta_cont_pre) :* xvar)
    M2_post 		= mean(w01 :* (y :- eta_cont_post) :* xvar)
    inf_cont_ps 	= lin_rep_ps * (M2_post :- M2_pre)'
    inf_cont 		= inf_cont :+ inf_cont_ps
    
	att_inf_func 	= ipw_att :+ inf_treat :- inf_cont

	st_store(.,rif,touse,att_inf_func)

	
 //  -15.80330618	
 //  9.087929526
 }
 
// Regression approach
  void reg_rc(string scalar y_, yy_, xvar_ , tmt_, trt_, wgt_ ,  touse, rif) {
    // main Loading variables
    real matrix y,  y00, y01,
				xvar, tmt, trt,   wgt
 	y    = st_data(.,y_, touse)
	yy   = st_data(.,yy_, touse)
	y00  = yy[,1]
	y01  = yy[,2]
		// verify xvar
 		if (xvar_==" ") {
			xvar=J(rows(y),1,1)	
		}
		else xvar=st_data(.,xvar_,touse),J(rows(y),1,1)
		
 	//xvar = st_data(.,xvar_, touse), J(rows(y),1,1)
	tmt  = st_data(.,tmt_, touse)
	trt  = st_data(.,trt_, touse)
	//psc  = logistic(st_data(.,pxb_, touse))
	//psv  = st_matrix(psv_)
	wgt  = st_data(.,wgt_, touse)
	
	real scalar nn
	nn = rows(y)
	
	real matrix w10, w11, w0
	
	w10 			= wgt :* trt :* (1 :- tmt)
    w11 			= wgt :* trt :* tmt
    w0 				= wgt :* trt
	
	w10				= w10:/mean(w10 )
	w11				= w11:/mean(w11 )
	w0				= w0:/mean(w0 )
	
	real matrix att_treat_pre, att_treat_post, att_cont,
				eta_treat_pre, eta_treat_post, eta_cont, reg_att
    att_treat_pre 	= w10 :* y
    att_treat_post 	= w11 :* y
    att_cont 		= w0 :* (y01 :- y00)
    eta_treat_pre 	= mean(att_treat_pre)
    eta_treat_post 	= mean(att_treat_post)
    eta_cont 		= mean(att_cont)
    reg_att 		= (eta_treat_post :- eta_treat_pre) :- eta_cont
	
	real matrix w_ols_pre, wols_eX_pre, XpX_inv_pre, lin_rep_ols_pre
    w_ols_pre 		= wgt :* (1 :- trt) :* (1 :- tmt)
    //wols_x_pre 	= w_ols_pre :* xvar
    wols_eX_pre 	= w_ols_pre :* (y :- y00) :* xvar
    XpX_inv_pre 	= invsym(quadcross(xvar,w_ols_pre, xvar)):*nn
    lin_rep_ols_pre = wols_eX_pre * XpX_inv_pre
	
	real matrix w_ols_post, wols_eX_post, XpX_inv_post, lin_rep_ols_post
    w_ols_post 		= wgt :* (1 :- trt) :* tmt
    //wols_x_post 	= w_ols_post :* xvar
    wols_eX_post 	= w_ols_post :* (y :- y01) :* xvar
    XpX_inv_post 	= invsym(quadcross(xvar, w_ols_post, xvar)):*nn
    lin_rep_ols_post = wols_eX_post * XpX_inv_post
    
	real matrix inf_treat, inf_cont_1, inf_cont_2_post, inf_cont_2_pre, inf_control
	//inf_treat_pre 	= (att_treat_pre :- w10 :* eta_treat_pre)
    //inf_treat_post 	= (att_treat_post :- w11 :* eta_treat_post)
    inf_treat 		= (att_treat_post :- w11 :* eta_treat_post) :- (att_treat_pre :- w10 :* eta_treat_pre)
    inf_cont_1 		= (att_cont :- w0 :* eta_cont)
    //M1 				= mean(w0 :* xvar)
    inf_cont_2_post = lin_rep_ols_post * mean(w0 :* xvar)'
    inf_cont_2_pre 	= lin_rep_ols_pre  * mean(w0 :* xvar)'
    inf_control 	= (inf_cont_1 :+ inf_cont_2_post :- inf_cont_2_pre)
    
	real matrix att_inf_func
	att_inf_func 	= reg_att :+ (inf_treat :- inf_control)

	st_store(.,rif,touse,att_inf_func)
	// Variance should be divided by n not n-1
	
 }

 /// Abadies IPW
 
  void ipw_abadie_rc(string scalar y_, xvar_ , tmt_, trt_, psv_, pxb_, wgt_ ,  touse, rif){
    // main Loading variables
    real matrix y, 	xvar, tmt, trt, psv, psc, wgt
				
		y    = st_data(.,y_, touse)
		// verify xvar
		if (xvar_==" ") {
			xvar=J(rows(y),1,1)	
		}
		else xvar=st_data(.,xvar_,touse),J(rows(y),1,1)			
		//xvar = st_data(.,xvar_, touse), J(rows(y),1,1)
		tmt  = st_data(.,tmt_, touse)
		trt  = st_data(.,trt_, touse)
		psc  = logistic(st_data(.,pxb_, touse))
		wgt  = st_data(.,wgt_, touse)
		psv  = st_matrix(psv_)
		real scalar nn
		nn = rows(y)
	
	real matrix w10, w11, w00, w01
    w10 				= wgt :* trt :* (1 :- tmt)
    w11 				= wgt :* trt :* tmt
    w00 				= wgt :* psc :* (1 :- trt) :* (1 :- tmt):/(1 :- psc)
    w01 				= wgt :* psc :* (1 :- trt) :* tmt:/(1 :- psc)
	
	real matrix pi_hat, lambda_hat, one_lambda_hat
	pi_hat       		= mean(wgt :* trt)
    lambda_hat 			= mean(wgt :* tmt)
    one_lambda_hat 		= mean(wgt :* (1 :- tmt))
	
	real matrix att_treat_pre, att_treat_post, att_cont_pre, att_cont_post,
				eta_treat_pre, eta_treat_post, eta_cont_pre, eta_cont_post
				
    att_treat_pre 		= w10 :* y:/(pi_hat :* one_lambda_hat)
    att_treat_post 		= w11 :* y:/(pi_hat :* lambda_hat)
    att_cont_pre 		= w00 :* y:/(pi_hat :* one_lambda_hat)
    att_cont_post 		= w01 :* y:/(pi_hat :* lambda_hat)
    eta_treat_pre 		= mean(att_treat_pre)
    eta_treat_post 		= mean(att_treat_post)
    eta_cont_pre 		= mean(att_cont_pre)
    eta_cont_post 		= mean(att_cont_post)
	
	real matrix ipw_att
    ipw_att 			= (eta_treat_post - eta_treat_pre) - (eta_cont_post - eta_cont_pre)
	
	real matrix lin_rep_ps, inf_treat_post, inf_treat_pret, inf_cont_post, inf_cont_pret
    //score_ps 			= wgt :* (trt :- psc) :* xvar
    //Hessian_ps 			= psv :* nn
    lin_rep_ps 			= (wgt :* (trt :- psc) :* xvar) * (psv :* nn)
	
	inf_treat_post 		=(att_treat_post :- eta_treat_post) :-
						((wgt :* trt :- pi_hat)     :* eta_treat_post:/pi_hat) :-
						((wgt :* tmt :- lambda_hat) :* eta_treat_post:/lambda_hat)
    ///inf_treat_post 		= inf_treat_post1 :+ inf_treat_post2 :+ inf_treat_post3
	
    inf_treat_pret 		= att_treat_pre :- eta_treat_pre :-
						 (wgt :*       trt  :- pi_hat) :* eta_treat_pre:/pi_hat :-
						 (wgt :* (1 :- tmt) :- one_lambda_hat) :* eta_treat_pre:/one_lambda_hat
    // inf_treat_pret 		= inf_treat_pre1 :+ inf_treat_pre2 :+ inf_treat_pre3
	
    inf_cont_post		= att_cont_post :- eta_cont_post :-
						  (wgt :* trt :- pi_hat) :* eta_cont_post:/pi_hat :-
						  (wgt :* tmt :- lambda_hat) :* eta_cont_post:/lambda_hat
    ///inf_cont_post = inf_cont_post1 :+ inf_cont_post2 :+ inf_cont_post3
    inf_cont_pret       = att_cont_pre :- eta_cont_pre :-
						  (wgt :* trt :- pi_hat) :* eta_cont_pre:/pi_hat :-
						  (wgt :* (1 :- tmt) :- one_lambda_hat) :* eta_cont_pre:/one_lambda_hat
     		
	//inf_cont_pret= inf_cont_pre1 :+ inf_cont_pre2 :+ inf_cont_pre3
	real matrix inf_logit, att_inf_func
    //mom_logit_pre 		= mean(-att_cont_pre :* xvar)
    //mom_logit_pre 		= mean(mom_logit_pre)
    //mom_logit_post 		= mean(-att_cont_post :* xvar)
    //mom_logit_post 		= mean(mom_logit_post)
    inf_logit 			= lin_rep_ps * (mean(-att_cont_post :* xvar) :-
										mean(-att_cont_pre :* xvar))'
	
    att_inf_func 		= ipw_att :+ (inf_treat_post :- inf_treat_pret) :- (inf_cont_post :- inf_cont_pret) :+ inf_logit
	
	st_store(.,rif,touse,att_inf_func)
		
//  -19.89330192
//  53.86822411
 }
 
/// DRDID_RC
 void drdid_rc(string scalar y_, yy_, xvar_ , tmt_, trt_, psv_, pxb_, wgt_ ,  touse, rif){
    // main Loading variables
    real matrix y,  y00, y01, y10, y11,
				xvar, tmt, trt, psv, psc, wgt
				
	y    = st_data(.,y_, touse)
	yy   = st_data(.,yy_, touse)
	y00  = yy[,1]
	y01  = yy[,2]
	y10  = yy[,3]
	y11  = yy[,4]
		// verify xvar
		if (xvar_==" ") {
			xvar=J(rows(y),1,1)	
		}
		else {
		    xvar=st_data(.,xvar_,touse),J(rows(y),1,1)	
		}
		
	//xvar = st_data(.,xvar_, touse), J(rows(y),1,1)
	tmt  = st_data(.,tmt_, touse)
	trt  = st_data(.,trt_, touse)
	psc  = logistic(st_data(.,pxb_, touse))
	wgt  = st_data(.,wgt_, touse)
	psv  = st_matrix(psv_)
	
	real matrix y0
	real scalar nn
	y0   = y00:*(-tmt:+1) + y01:*tmt
	nn = rows(y)
	
	real matrix w10, w11, w00, w01, w1
    w10 = wgt :* trt :* (1 :- tmt)
    w11 = wgt :* trt :* tmt
    w00 = wgt :* psc :* (1 :- trt) :* (1 :- tmt):/(1 :- psc)
    w01 = wgt :* psc :* (1 :- trt) :* tmt:/(1 :- psc)
	w1  = wgt :* trt
	
	w00 = w00:/mean(w00 )
	w01 = w01:/mean(w01 )
	w10 = w10:/mean(w10 )
	w11 = w11:/mean(w11 )
	w1 = w1:/mean(w1 )
	
	real matrix att_treat_pre, att_treat_post,  att_cont_pre, att_cont_post,
				att_trt_post, att_trtt1_post, att_trt_pre, att_trtt0_pre,
				eta_treat_pre, eta_treat_post,  eta_cont_pre, eta_cont_post,
				eta_trt_post, eta_trtt1_post, eta_trt_pre, eta_trtt0_pre
				
    att_treat_pre 		= w10 :* (y :- y0)
    att_treat_post 		= w11 :* (y :- y0)
    att_cont_pre  		= w00 :* (y :- y0)
    att_cont_post  		= w01 :* (y :- y0)
	
    att_trt_post   		= w1  :* (y11 :- y01)
    att_trtt1_post 		= w11 :* (y11 :- y01)
    att_trt_pre   		= w1  :* (y10 :- y00)
    att_trtt0_pre 		= w10 :* (y10 :- y00)
	
    eta_treat_pre 		= mean(att_treat_pre)
    eta_treat_post 		= mean(att_treat_post)
    eta_cont_pre  		= mean(att_cont_pre)
    eta_cont_post  		= mean(att_cont_post)
    eta_trt_post   		= mean(att_trt_post)
    eta_trtt1_post 		= mean(att_trtt1_post)
    eta_trt_pre   		= mean(att_trt_pre)
    eta_trtt0_pre 		= mean(att_trtt0_pre)
	
	real matrix trtr_att
    trtr_att      		= (eta_treat_post :- eta_treat_pre) :- (eta_cont_post :- eta_cont_pre) :+ (eta_trt_post :- eta_trtt1_post) :- (eta_trt_pre :- eta_trtt0_pre)
	
	real matrix wgt_ols_pre, XpX_inv_pre, lin_ols_pre, 
				wgt_ols_post, XpX_inv_post, lin_ols_post,
				XpX_inv_pre_treat, lin_ols_pre_treat,
				XpX_inv_post_treat, lin_ols_post_treat
	wgt_ols_pre     	= wgt :* (1 :- trt) :* (1 :- tmt)
    //wols_x_pre          = wgt_ols_pre :* xvar
    //wols_eX_pre         = wgt_ols_pre :* (y :- y00) :* xvar
    XpX_inv_pre         = invsym(quadcross(xvar,wgt_ols_pre, xvar)):*nn
    lin_ols_pre 		= ( wgt_ols_pre :* (y :- y00) :* xvar) * XpX_inv_pre
    
	wgt_ols_post     	= wgt :* (1 :- trt) :* tmt
    //wols_x_post         = wgt_ols_post :* xvar
    //wols_eX_post        = wgt_ols_post :* (y :- y01) :* xvar
    XpX_inv_post 		= invsym(quadcross(xvar,wgt_ols_post, xvar)):*nn
	lin_ols_post 		= (wgt_ols_post :* (y :- y01) :* xvar) * XpX_inv_post
	    
    //wols_x_pre_treat 	= w10 :* xvar
    //wols_eX_pre_treat 	= w10 :* (y :- y10) :* xvar
    XpX_inv_pre_treat 	= invsym(quadcross(xvar, w10, xvar)):*nn
	lin_ols_pre_treat 	= ( w10 :* (y :- y10) :* xvar) * XpX_inv_pre_treat
	
	//wols_x_post_treat   = w11 :* xvar
    //wols_eX_post_treat  = w11 :* (y :- y11) :* xvar
    XpX_inv_post_treat  = invsym(quadcross(xvar, w11, xvar)):*nn
    lin_ols_post_treat 	= (w11 :* (y :- y11) :* xvar) * XpX_inv_post_treat
	
	real matrix lin_rep_ps, inf_treat_pre, inf_treat_post
    // check psv for probit
	//score_ps 			= wgt :* (trt :- psc) :* xvar
    //Hessian_ps 			= psv :* nn
    lin_rep_ps 			= (wgt :* (trt :- psc) :* xvar) * (psv :* nn)
    inf_treat_pre 		= att_treat_pre  :- w10 :* eta_treat_pre 
    inf_treat_post 		= att_treat_post :- w11 :* eta_treat_post
	
	real matrix M1_post, M1_pre, inf_treat_or_post, inf_treat_or_pre
	
    M1_post 			= -mean(w11 :* tmt :* xvar)
    M1_pre 				= -mean(w10 :* (1 :- tmt) :* xvar)
	inf_treat_or_post 	= lin_ols_post * M1_post'
    inf_treat_or_pre 	= lin_ols_pre * M1_pre'
	
	real matrix inf_treat_or, inf_treat, inf_cont_pre, inf_cont_post
	
    inf_treat_or 		= inf_treat_or_post :+ inf_treat_or_pre
    inf_treat 			= inf_treat_post :- inf_treat_pre :+ inf_treat_or
    inf_cont_pre 		= att_cont_pre  :- w00 :* eta_cont_pre 
    inf_cont_post 		= att_cont_post :- w01 :* eta_cont_post 
    
	real matrix M2_pre, M2_post, inf_cont_ps, M3_post, M3_pre, inf_cont_or_post, inf_cont_or_pre
	M2_pre 				= mean(w00 :* (y :- y0 :- eta_cont_pre) :* xvar)
    M2_post 			= mean(w01 :* (y :- y0 :- eta_cont_post) :* xvar)
    inf_cont_ps 		= lin_rep_ps * (M2_post :- M2_pre)'
	
    M3_post 			= -mean(w01 :* tmt :* xvar)
    M3_pre 				= -mean(w00 :* (1 :- tmt) :* xvar)
    inf_cont_or_post 	= lin_ols_post * M3_post'
    inf_cont_or_pre 	= lin_ols_pre  * M3_pre'
	
	real matrix inf_cont_or, inf_cont, trtr_eta_inf_func1
	
    inf_cont_or 		= inf_cont_or_post :+ inf_cont_or_pre
    inf_cont 			= inf_cont_post    :- inf_cont_pre :+ inf_cont_ps :+ inf_cont_or
    trtr_eta_inf_func1 	= inf_treat :- inf_cont
	
	real matrix inf_eff, mom_post, mom_pre, inf_or
    //inf_eff1 			= att_trt_post   :- w1  :* eta_trt_post   
    //inf_eff2 			= att_trtt1_post :- w11 :* eta_trtt1_post 
    //inf_eff3 			= att_trt_pre    :- w1  :* eta_trt_pre   
    //inf_eff4 			= att_trtt0_pre  :- w10 :* eta_trtt0_pre 
    inf_eff 			= ((att_trt_post   :- w1  :* eta_trt_post)    :- 
						   (att_trtt1_post :- w11 :* eta_trtt1_post)) :-
						  ((att_trt_pre    :- w1  :* eta_trt_pre)     :- 
						   (att_trtt0_pre :- w10 :* eta_trtt0_pre))
    mom_post 			= mean((w1 :- w11) :* xvar)
    mom_pre 			= mean((w1 :- w10) :* xvar)
	// check this
    //inf_or_post 		= ((lin_ols_post_treat :- lin_ols_post) * mom_post')
    //inf_or_pre 		= 
    inf_or 				= ((lin_ols_post_treat :- lin_ols_post) * mom_post') :- ((lin_ols_pre_treat :- lin_ols_pre) * mom_pre')
	
	att_inf_func	 	= trtr_att :+ trtr_eta_inf_func1 :+ inf_eff :+ inf_or
	
	st_store(.,rif,touse,att_inf_func)
	// Variance should be divided by n not n-1
 	
//  -.1677954483	
// .2008991705	
 }

/// DRDID_RC1
void drdid_rc1(string scalar y_, yy_, xvar_ , tmt_, trt_, psv_, pxb_, wgt_ ,  touse, rif){
    // main Loading variables
    real matrix y,  y00, y01, y10, y11,
				xvar, tmt, trt, psv, psc, wgt
				
	y    = st_data(.,y_, touse)
	yy   = st_data(.,yy_, touse)
	y00  = yy[,1]
	y01  = yy[,2]
	y10  = yy[,3]
	y11  = yy[,4]
		// verify xvar
		if (xvar_==" ") {
			xvar=J(rows(y),1,1)	
		}
		else xvar=st_data(.,xvar_,touse),J(rows(y),1,1)		
	//xvar = st_data(.,xvar_, touse), J(rows(y),1,1)
	
	tmt  = st_data(.,tmt_, touse)
	trt  = st_data(.,trt_, touse)
	psc  = logistic(st_data(.,pxb_, touse))
	wgt  = st_data(.,wgt_, touse)
	psv  = st_matrix(psv_)
	
	real matrix y0
	real scalar nn
	y0   = y00:*(-tmt:+1) + y01:*tmt
	nn = rows(y)
    
	real matrix w10, w11, w00, w01
    w10 		= wgt :* trt :* (1 :- tmt)
    w11 		= wgt :* trt :* tmt
    w00 		= wgt :* psc :* (1 :- trt) :* (1 :- tmt):/(1 :- psc)
    w01 		= wgt :* psc :* (1 :- trt) :* tmt:/(1 :- psc)
	
	w00 = w00:/mean(w00 )
	w01 = w01:/mean(w01 )
	w10 = w10:/mean(w10 )
	w11 = w11:/mean(w11 )
	
	real matrix eta_trt_pre, eta_trt_post, eta_cont_pre, eta_cont_post,
				att_trt_pre, att_trt_post, att_cont_pre, att_cont_post
				
    eta_trt_pre 	= w10 :* (y :- y0)
    eta_trt_post 	= w11 :* (y :- y0)
    eta_cont_pre 	= w00 :* (y :- y0)
    eta_cont_post 	= w01 :* (y :- y0)
    att_trt_pre 	= mean(eta_trt_pre)
    att_trt_post 	= mean(eta_trt_post)
    att_cont_pre 	= mean(eta_cont_pre)
    att_cont_post 	= mean(eta_cont_post)
	
	real matrix trtr_att, w_ols_pre, wols_eX_pre, lin_rep_ols_pre,
						  w_ols_post, wols_eX_post, lin_rep_ols_post
    trtr_att 		= (att_trt_post :- att_trt_pre) :- (att_cont_post :- att_cont_pre)
	
    w_ols_pre 		= wgt :* (1 :- trt) :* (1 :- tmt)
    wols_eX_pre 	= w_ols_pre :* (y :- y00) :* xvar
 
    lin_rep_ols_pre = wols_eX_pre * invsym(quadcross(xvar, w_ols_pre, xvar)):*nn
	
    w_ols_post  	= wgt :* (1 :- trt) :* tmt
    wols_eX_post 	= w_ols_post  :* (y :- y01) :* xvar
    lin_rep_ols_post = wols_eX_post * invsym(quadcross( xvar,w_ols_post, xvar)):*nn
	
	real matrix lin_rep_ps, inf_trt_pre, inf_trt_post, M1_post, M1_pre
    //score_ps 		= wgt :* (trt :- psc) :* xvar
    //Hessian_ps 		= psv :* nn
    lin_rep_ps 		= (wgt :* (trt :- psc) :* xvar) * (psv :* nn)
    inf_trt_pre 	= eta_trt_pre :- w10 :* att_trt_pre
    inf_trt_post 	= eta_trt_post :- w11 :* att_trt_post
	M1_post 		= -mean(w11 :* tmt :* xvar)
    M1_pre 			= -mean(w10 :* (1 :- tmt) :* xvar)
    
	real matrix inf_trt_or, inf_trt
    inf_trt_or 		= (lin_rep_ols_post * M1_post') :+ (lin_rep_ols_pre * M1_pre')
    inf_trt 		= inf_trt_post :- inf_trt_pre :+ inf_trt_or
    
	real matrix inf_cont_pre , inf_cont_post, M2_pre,  M2_post, inf_cont_ps, M3_post, M3_pre

	inf_cont_pre 	= eta_cont_pre :- w00 :* att_cont_pre
    inf_cont_post 	= eta_cont_post :- w01 :* att_cont_post
    M2_pre 			= mean(w00 :* (y :- y0 :- att_cont_pre) :* xvar)
    M2_post 		= mean(w01 :* (y :- y0 :- att_cont_post) :* xvar)
    inf_cont_ps 	= lin_rep_ps * (M2_post :- M2_pre)'
    M3_post 		= -mean(w01 :* tmt :* xvar)
    M3_pre 			= -mean(w00 :* (1 :- tmt) :* xvar)

	real matrix inf_cont_or, inf_cont, att_inf_func
    inf_cont_or 	= (lin_rep_ols_post * M3_post') :+ (lin_rep_ols_pre * M3_pre')
    inf_cont 		= inf_cont_post :- inf_cont_pre :+ inf_cont_ps :+ inf_cont_or
    att_inf_func 	= trtr_att :+ inf_trt :- inf_cont
	
	st_store(.,rif,touse,att_inf_func)
	// Variance should be divided by n not n-1
 	
 // -3.633433441	
//3.107123089
}
 
/// drdid_imp_rc
void drdid_imp_rc(string scalar y_, yy_, xvar_ , tmt_, trt_, psv_, pxb_, wgt_ ,  touse, rif){
    // main Loading variables
    real matrix y,  y00, y01, y10, y11,
				xvar, tmt, trt, psv, psc, wgt
				
	y    = st_data(.,y_, touse)
	yy   = st_data(.,yy_, touse)
	y00  = yy[,1]
	y01  = yy[,2]
	y10  = yy[,3]
	y11  = yy[,4]
		// verify xvar
		if (xvar_==" ") {
			xvar=J(rows(y),1,1)	
		}
		else xvar=st_data(.,xvar_,touse),J(rows(y),1,1)		
//	xvar = st_data(.,xvar_, touse), J(rows(y),1,1)
	tmt  = st_data(.,tmt_, touse)
	trt  = st_data(.,trt_, touse)
	psc  = logistic(st_data(.,pxb_, touse))
	wgt  = st_data(.,wgt_, touse)
	psv  = st_matrix(psv_)
	
	real matrix y0
	real scalar nn
	y0   = y00:*(-tmt:+1) + y01:*tmt
	nn = rows(y)

	real matrix w10, w11, w00, w01, w1
    w10 = wgt :* trt :* (1 :- tmt)
    w11 = wgt :* trt :* tmt
    w00 = wgt :* psc :* (1 :- trt) :* (1 :- tmt):/(1 :-  psc)
    w01 = wgt :* psc :* (1 :- trt) :* tmt:/(1 :- psc)
    w1  = wgt :* trt

	w00 = w00:/mean(w00 )
	w01 = w01:/mean(w01 )
	w10 = w10:/mean(w10 )
	w11 = w11:/mean(w11 )
	w1 = w1:/mean(w1 )
	
	real matrix att_treat_pre, att_treat_post, att_cont_pre, att_cont_post, att_trt_post, att_trtt1_post,
				att_trt_pre, att_trtt0_pre
	
    att_treat_pre  = w10 :* (y :- y0)
    att_treat_post = w11 :* (y :- y0)
    att_cont_pre   = w00 :* (y :- y0)
    att_cont_post  = w01 :* (y :- y0)
    att_trt_post   = w1  :* (y11 :- y01)
    att_trtt1_post = w11 :* (y11 :- y01)
    att_trt_pre    = w1  :* (y10 :- y00)
    att_trtt0_pre  = w10 :* (y10 :- y00)
	
	real matrix eta_treat_pre, eta_treat_post, eta_cont_pre, eta_cont_post, eta_trt_post, eta_trtt1_post,
				eta_trt_pre, eta_trtt0_pre
				
    eta_treat_pre  = mean(att_treat_pre)
    eta_treat_post = mean(att_treat_post)
    eta_cont_pre   = mean(att_cont_pre)
    eta_cont_post  = mean(att_cont_post)
    eta_trt_post   = mean(att_trt_post)
    eta_trtt1_post = mean(att_trtt1_post)
    eta_trt_pre    = mean(att_trt_pre)
    eta_trtt0_pre  = mean(att_trtt0_pre)
	
	real matrix trtr_att
	
    trtr_att       = (eta_treat_post :- eta_treat_pre) :- (eta_cont_post :- eta_cont_pre) :+ (eta_trt_post :- eta_trtt1_post) :- (eta_trt_pre :- eta_trtt0_pre)
	
	real matrix inf_treat,  inf_cont,  att_inf_func1,  inf_eff, att_inf_func

	
	inf_treat      = (att_treat_post :- w11 :* eta_treat_post) :- (att_treat_pre  :- w10 :* eta_treat_pre)
    inf_cont       = (att_cont_post  :- w01 :* eta_cont_post)  :- (att_cont_pre   :- w00 :* eta_cont_pre)
	att_inf_func1  = inf_treat :- inf_cont
    inf_eff        =  ((att_trt_post   :- w1  :* eta_trt_post) :- 
					   (att_trtt1_post :- w11 :* eta_trtt1_post)) :- 
					  ((att_trt_pre    :- w1  :* eta_trt_pre) :- 
					   (att_trtt0_pre  :- w10 :* eta_trtt0_pre))
	
    att_inf_func  = trtr_att :+ att_inf_func1 :+ inf_eff
	
	st_store(.,rif,touse,att_inf_func)
	// Variance should be divided by n not n-1
 	
//-.2088586355	
//.2003375215	
}

/// drdid_imp_rc1
void drdid_imp_rc1(string scalar y_, yy_, xvar_ , tmt_, trt_, psv_, pxb_, wgt_ ,  touse, rif){
    // main Loading variables
    real matrix y,  y00, y01, y10, y11,
				xvar, tmt, trt, psv, psc, wgt
	y    = st_data(.,y_, touse)
	yy   = st_data(.,yy_, touse)
	y00  = yy[,1]
	y01  = yy[,2]
	y10  = yy[,3]
	y11  = yy[,4]
		// verify xvar
		if (xvar_==" ") {
			xvar=J(rows(y),1,1)	
		}
		else xvar=st_data(.,xvar_,touse),J(rows(y),1,1)		
	//xvar = st_data(.,xvar_, touse), J(rows(y),1,1)
	tmt  = st_data(.,tmt_, touse)
	trt  = st_data(.,trt_, touse)
	psc  = logistic(st_data(.,pxb_, touse))
	wgt  = st_data(.,wgt_, touse)
	psv  = st_matrix(psv_)
	
	real matrix y0
	real scalar nn
	y0   = y00:*(-tmt:+1) + y01:*tmt
	nn = rows(y)
	
	real matrix w10, w11, w00, w01, 
				eta_treat_pre, eta_treat_post, eta_cont_pre, eta_cont_post,
				att_treat_pre, att_treat_post, att_cont_pre, att_cont_post
				
    w10 			= wgt :* trt :* (1 :- tmt)
    w11 			= wgt :* trt :* tmt
    w00 			= wgt :* psc:* (1 :- trt) :* (1 :- tmt):/(1 :-  psc)
	w01 			= wgt :* psc:* (1 :- trt) :* tmt:/(1 :- psc)
    
	eta_treat_pre 	= w10 :* (y :- y0):/mean(w10)
    eta_treat_post 	= w11 :* (y :- y0):/mean(w11)
    eta_cont_pre 	= w00 :* (y :- y0):/mean(w00)
    eta_cont_post 	= w01 :* (y :- y0):/mean(w01)
    
	att_treat_pre 	= mean(eta_treat_pre)
    att_treat_post 	= mean(eta_treat_post)
    att_cont_pre 	= mean(eta_cont_pre)
    att_cont_post 	= mean(eta_cont_post)
    
	real matrix trtr_att, inf_treat, inf_cont, att_inf_func
	trtr_att 		= (att_treat_post :- att_treat_pre) :- (att_cont_post :-  att_cont_pre)
    inf_treat 		= (eta_treat_post :- w11 :* att_treat_post:/mean(w11) ):- (eta_treat_pre :- w10 :* att_treat_pre:/mean(w10))
    inf_cont 		= (eta_cont_post :- w01 :* att_cont_post:/mean(w01)) :- ( eta_cont_pre :- w00 :* att_cont_pre:/mean(w00))
    att_inf_func 	= trtr_att :+ inf_treat :- inf_cont
	
	// Wrapping up
	st_store(.,rif,touse,att_inf_func)
	// Variance should be divided by n not n-1
 	
//  -3.683728719	
//	 3.114495585
}

// Clustered Standard errors

void clusterse(real matrix iiff, cl, V, real scalar cln){
    /// estimates Clustered Standard errors
    real matrix ord, xcros, ifp, info, vv 

	ord  = order(cl,1)
	iiff = iiff[ord,]
	cl   = cl[ord,]
	// check how I cleaned data!
	info  = panelsetup(cl,1)
	// faster Cluster? Need to do this for mmqreg
	ifp   = panelsum(iiff,info)
	xcros = quadcross(ifp,ifp)
	real scalar nt, nc
	nt=rows(iiff)
	nc=rows(info)
	V =	xcros/(nt^2)
	cln=nc
	// Esto es para ver como hacer clusters.
	//*nc/(nc-1)
	//st_matrix(V,    vv)
	//st_numscalar(ncl, nc)
	//        ^     ^
	//        |     |
	//      stata   mata
}
//** Simple Bootstrap.
 
real matrix mboot_did(real matrix rif, mean_rif, real scalar reps, wbtype) {
	real matrix yy, bsmean
	yy=rif:-mean_rif
 	bsmean=J(reps,cols(yy),0)
	real scalar i,n, k1, k2
	n=rows(yy)
	k1=((1+sqrt(5))/(2*sqrt(5)))
	k2=0.5*(1+sqrt(5)) 
	// WBootstrap:Mammen 
	if (wbtype==1) {			
		for(i=1;i<=reps;i++){
			bsmean[i,]=mean(yy:*(k2:-sqrt(5)*(rbinomial(n,1,1,k1))) )	
		}
	}
	
	else if (wbtype==2) {
		for(i=1;i<=reps;i++){
			bsmean[i,]=mean(yy:*(1:-2*rbinomial(n,1,1,0.5) ) )	
		}
	}
	
	return(bsmean)
}

real matrix mboot_didc(real matrix rif, mean_rif, clv , real scalar reps, wbtype, nc) {
	real matrix yy, bsmean
	real matrix sclv, wmult
	real scalar i,n, k1, k2, nn
	yy=rif:-mean_rif
 	bsmean=J(reps,cols(yy),0)
	
	nn=rows(uniqrows(clv))
	nc=nn
	k1=((1+sqrt(5))/(2*sqrt(5)))
	k2=0.5*(1+sqrt(5))
	if (wbtype==1) {		
 
		for(i=1;i<=reps;i++){
		    wmult=rbinomial(nn,1,1,k1)
			bsmean[i,]=mean(yy:*(k2:-sqrt(5)*wmult[clv] ) )	
			
		}
	}
	
	else if (wbtype==2) {
		for(i=1;i<=reps;i++){
		    wmult=(rbinomial(nn,1,1,0.5))
			bsmean[i,]=mean(yy:*(1:-2* wmult[clv] ) )	
		}
	}
	return(bsmean)
	
}

///clust("`att'","`cluster'","`touse'", "`V'",
///mboot("`att'", "`touse'", "`V'","`ciband'","`cluster'", `reps', `wbtype', `ci')
void make_tbl(string scalar iiff, clv, touse, // RIF, Cluster variable, and sample
			  bb_ , VV_, ccbb, cln,                   // WHere to save bb and VV. Also #clv
			  wboot,                          // options for Wbootstrap. Reps, wbtype and CI
			  real scalar reps, wbtype, ci ){
			      
	real matrix rif, clvar , bb, VV, ord, info, ifp
	real matrix cband
	real scalar nobs, nc
	
	// nc number of clusters
	rif = st_data(.,iiff,touse)
	  
	// real scalar cln
	bb  = mean(rif)
	nobs     = rows(rif)
	// simple
 
	if ((clv==" ") & (wboot==" ")) {	
		VV=quadcrossdev(rif,bb,rif,bb):/ (nobs^2) 
	}
	
	// cluster std
	if ((clv!=" ") & (wboot==" ")) {
		clvar=st_data(.,clv,touse)
		ord   = order(clvar,1)
		rif   = rif[ord,]
		clvar = clvar[ord,]
		
		info  = panelsetup(clvar,1)
 		ifp   = panelsum((rif:-bb),info)
 		VV    = quadcross(ifp,ifp)/nobs:^2
		nc=rows(info)		
	}
	
	// wboot no cluster
	
	if (wboot!=" ") {
	    
		mboot(rif,bb, VV, cband, clv, touse ,reps, wbtype, ci, nc )
	}
	
	st_matrix(bb_,bb)
	st_matrix(VV_,VV)
	st_matrix(ccbb,cband)
	st_numscalar(cln, nc)
	 
	
 } 
 
real matrix iqrse(real matrix y) {
    real scalar q25,q75
	q25=ceil((rows(y)+1)*.25)
	q75=ceil((rows(y)+1)*.75)
	real scalar j
	real matrix iqrs
	iqrs=J(1,cols(y),0)
	for(j=1;j<=cols(y);j++){
	    y=sort(y,j)
		iqrs[,j]=(y[q75,j]-y[q25,j]):/(invnormal(.75)-invnormal(.25) )
	}
	return(iqrs)
}

real vector qtp(real matrix y, real scalar p) {
    real scalar k, i, q
	real matrix yy, qq
	qq=J(1,0,.)
	k = cols(y)
	for(i=1;i<=k;i++){
		yy=sort(y[,i],1)
		q=ceil((rows(yy)+1)*p)    
		qq=qq,yy[q,]
	}
    
	return(qq)
}

void mboot(real matrix rif,mean_rif, vv, cband, string scalar clv, touse,  real scalar reps, wbtype, ci , nc ) {
    //, real scalar reps, wbtype, ci 
    real matrix fr, qt
	//real scalar reps, wbtype
 
	real matrix ifse  
	// this gets the Bootstraped values
	if (clv ==" ") {
		fr=mboot_did(rif,mean_rif      , reps, wbtype)
		ifse = iqrse(fr)
		qt = qtp(abs(fr :/ ifse),ci)
		cband=( mean_rif':-qt':* ifse' ,  
				mean_rif':+qt':* ifse' , mean_rif, ifse, qt' )
	}
	else {
 		clvar=st_data(.,clv, touse)
		fr=mboot_didc(rif,mean_rif, clvar, reps, wbtype, nc )
		ifse = iqrse(fr)
		qt = qtp(abs(fr :/ ifse),ci)
		cband=( mean_rif':-qt':* ifse' ,  
				mean_rif':+qt':* ifse' , mean_rif, ifse, qt' )
	}
	vv=quadcross(ifse,ifse):*I(cols(ifse))

}
 
 
/// qtp(abs(xx/ iqrse(xx)),.95) 
end

program addr, rclass
	return `0'
end 
